This notebook is meant for release to the general public, to reproduce the figures from our paper "Multiple mechanisms mediate the suppression of motion vision during escape maneuvers in flying Drosophila" https://doi.org/10.1101/2022.02.03.478949
These imports are necessary to run the notebook. The production environment is listed in "working_environment.txt".
The custom pathdefs package is a convenience for the author to simplify running on different machines. It exports two variables SOURCEPATH and DATAPATH , which may just as well be declared as regular variables.
The paths where the respective files are to be found will be different for other users as well, of course. Adapt to your folder structure.
%matplotlib notebook
import os
from pathdefs import *
os.chdir(SOURCEPATH)
import importlib
import Looming_with_background_Analyzer as lwb
#import baseClasses as bc
from spontaneousAnalyzer import DummyAnalyzer
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy import stats
from scipy import ndimage
import glob
import helperLibrary
# import re
# import copy
optomotorColor = np.array([75,75,75])/255
optomotorAlpha = 0.2
plotColors = [(.0, 100/255, 178/255), # blueish
(54/255, 170/255, 80/255), # some green
(200/255, 74/255, .0), # orange-red
optomotorColor]
plotAlphas = [0.3, # same order as colors
0.38,
0.31,
optomotorAlpha]
PLOT_TIMEWINDOW = [-0.65,0.46] # for most plots
# DEPENDENT ON PRECISE FOLDER STRUCTURE, PROBABLY DIFFERENT FOR ANYONE ELSE
# ADJUST THIS!!
saccLeft_HyPol = [('vonBettina/25072018_CS', 'Fly1_WTCS_age_3d_cell1_'),
('vonBettina/13022019_G06', 'Fly1_U29xG06_age_3d_cell1_'),
('vonBettina/26072019_CS', 'Fly1_WTCS_age_3d_cell1_'),
('2020_01_29', 'Fly3_U29xG16_age_5d_cell1_'),
('2020_03_16', 'Fly1_X20a_age_3d_cell1_'),
('2020_04_20', 'fly4_X21_age_3d_cell1_') ]
saccLeft_DePol = [('vonBettina/06112019_CS', 'Fly1_WTCS_age_3d_cell1_'),
('2020_01_29', 'Fly1_U29xG16_age_5d_cell1_') ]
saccLeft_NoPol = [('vonBettina/04112019_G06', 'Fly1_G06xU29_age_3d_cell1_'),
('2019_12_09', 'Fly1_U29xG16_age_2d_cell2_'),
('2020_04_15', 'Fly3_X21_age_3d_cell1_'),
('2020_04_17', 'Fly3_X21_age_5d_cell1_'),
('2021_02_17', 'fly2_X20_age_2d_cell1_') ]
saccRight_HyPol = [('vonBettina/25072018_CS', 'Fly1_WTCS_age_3d_cell1_'),
('vonBettina/13022019_G06', 'Fly1_U29xG06_age_3d_cell1_'),
('vonBettina/26072019_CS', 'Fly1_WTCS_age_3d_cell1_')]
saccRight_DePol = [('2020_03_16', 'Fly1_X20a_age_3d_cell1_'),
('vonBettina/04112019_G06', 'Fly1_G06xU29_age_3d_cell1_'),
('2020_01_29', 'Fly1_U29xG16_age_5d_cell1_'),
('2020_01_29', 'Fly3_U29xG16_age_5d_cell1_') ]
saccRight_NoPol = [('vonBettina/06112019_CS', 'Fly1_WTCS_age_3d_cell1_'), #(?)
('2019_12_09', 'Fly1_U29xG16_age_2d_cell2_'),
('2020_04_15', 'Fly3_X21_age_3d_cell1_'),
('2020_04_17', 'Fly3_X21_age_5d_cell1_'),
('2020_04_20', 'fly4_X21_age_3d_cell1_'),
('2021_02_17', 'fly2_X20_age_2d_cell1_') ]
# refer to https://gist.github.com/dmeliza/3251476
# dmeliza dmeliza/scalebars.py
from matplotlib.offsetbox import AnchoredOffsetbox
class AnchoredScaleBar(AnchoredOffsetbox):
def __init__(self, transform, sizex=0, sizey=0, labelx=None, labely=None, loc=4,
pad=0.1, borderpad=0.1, sep=2, prop=None, barcolor="black", barwidth=None,
**kwargs):
"""
Draw a horizontal and/or vertical bar with the size in data coordinate
of the give axes. A label will be drawn underneath (center-aligned).
- transform : the coordinate frame (typically axes.transData)
- sizex,sizey : width of x,y bar, in data units. 0 to omit
- labelx,labely : labels for x,y bars; None to omit
- loc : position in containing axes
- pad, borderpad : padding, in fraction of the legend font size (or prop)
- sep : separation between labels and bars in points.
- **kwargs : additional arguments passed to base class constructor
"""
from matplotlib.patches import Rectangle
from matplotlib.offsetbox import AuxTransformBox, VPacker, HPacker, TextArea, DrawingArea
bars = AuxTransformBox(transform)
if sizex:
bars.add_artist(Rectangle((0,0), sizex, 0, ec=barcolor, lw=barwidth, fc="none"))
if sizey:
bars.add_artist(Rectangle((0,0), 0, sizey, ec=barcolor, lw=barwidth, fc="none"))
if sizex and labelx:
self.xlabel = TextArea(labelx, minimumdescent=False, textprops={'fontname':'arial', 'fontsize':8.})
bars = VPacker(children=[bars, self.xlabel], align="center", pad=0, sep=sep)
if sizey and labely:
self.ylabel = TextArea(labely, textprops={'fontname':'arial', 'fontsize':8.})
bars = HPacker(children=[self.ylabel, bars], align="center", pad=0, sep=sep)
AnchoredOffsetbox.__init__(self, loc, pad=pad, borderpad=borderpad,
child=bars, prop=prop, frameon=False, **kwargs)
def add_scalebar(ax, matchx=True, matchy=True, hidex=True, hidey=True, **kwargs):
""" Add scalebars to axes
Adds a set of scale bars to *ax*, matching the size to the ticks of the plot
and optionally hiding the x and y axes
- ax : the axis to attach ticks to
- matchx,matchy : if True, set size of scale bars to spacing between ticks
if False, size should be set using sizex and sizey params
- hidex,hidey : if True, hide x-axis and y-axis of parent
- **kwargs : additional arguments passed to AnchoredScaleBars
Returns created scalebar object
"""
def f(axis):
l = axis.get_majorticklocs()
return len(l)>1 and (l[1] - l[0])
if matchx:
kwargs['sizex'] = f(ax.xaxis)
kwargs['labelx'] = str(kwargs['sizex'])
if matchy:
kwargs['sizey'] = f(ax.yaxis)
kwargs['labely'] = str(kwargs['sizey'])
sb = AnchoredScaleBar(ax.transData, **kwargs)
ax.add_artist(sb)
if hidex : ax.xaxis.set_visible(False)
if hidey : ax.yaxis.set_visible(False)
if hidex and hidey: ax.set_frame_on(False)
return sb
saSacc_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_straight.neomat'))
# keep only one fly
saSacc_Rhyp.clean_flyIDs()
rm = np.where(saSacc_Rhyp.fly_indices!=3)[0]
saSacc_Rhyp.remove_duds_byList(rm)
saStrg_Rhyp.clean_flyIDs()
rm = np.where(saStrg_Rhyp.fly_indices!=3)[0]
saStrg_Rhyp.remove_duds_byList(rm)
saSacc_Rhyp.clean_flyIDs()
saStrg_Rhyp.clean_flyIDs()
saSacc_Rhyp.reallocate_fly_indices()
saStrg_Rhyp.reallocate_fly_indices()
print(saSacc_Rhyp.flyID)
print(saStrg_Rhyp.flyID)
{0: 'fly1_WTCS_20190726_cell1'}
{0: 'fly1_WTCS_20190726_cell1'}
# keep only looming left trials
stimID = 2
saStrg_Rhyp.remove_duds_byList(np.where(saStrg_Rhyp.group_indices!=stimID)[0])
saSacc_Rhyp.remove_duds_byList(np.where(saSacc_Rhyp.group_indices!=stimID)[0])
fig, axes = plt.subplots(3,1, figsize=(5,8),sharex='all', gridspec_kw={'height_ratios': [2.5,2.5,1]})
saSacc = saSacc_Rhyp
saStrg = saStrg_Rhyp
f,ax = fig, axes[1]
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom+PLOT_TIMEWINDOW[0]*saSacc.samplingRate), int(maxLoom+PLOT_TIMEWINDOW[1]*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate+PLOT_TIMEWINDOW[0]
d = saStrg.segments[0].data[saStrg.channelmap['ephys']]*100. # factor of 100 converts to mV
ax.plot(timeSlice, d[timeLims[0]:timeLims[1]], color=plotColors[0])
d = saSacc.segments[1].data[saSacc.channelmap['ephys']]*100.
ax.plot(timeSlice, d[timeLims[0]:timeLims[1]], color=plotColors[2])
ax.axvline(0, linestyle='--', color='k')
ax.axhline(0, linestyle=':', color='k')
# ax.legend(['no saccade','with saccade'], font='Arial')
ax.set_ylabel('U [mV]', fontname='Arial')
plt.tight_layout()
ax.xaxis.set_visible(False)
ax.set_xlim(PLOT_TIMEWINDOW)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['right'].set_visible(False)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, hidey=False, sizey=0.0, barwidth=3)
plt.tight_layout()
saSacc = saSacc_Rhyp
saStrg = saStrg_Rhyp
f,ax = fig, axes[0]
# f,ax=plt.subplots(1,1)
preSacc = 0.8 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom+PLOT_TIMEWINDOW[0]*saSacc.samplingRate), int(maxLoom+PLOT_TIMEWINDOW[1]*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate+PLOT_TIMEWINDOW[0]
d = saStrg.segments[0].data[saStrg.channelmap['wingDiff']]*100.
ax.plot(timeSlice, d[timeLims[0]:timeLims[1]], color=plotColors[0])
d = saSacc.segments[1].data[saSacc.channelmap['wingDiff']]*100.
ax.plot(timeSlice, d[timeLims[0]:timeLims[1]], color=plotColors[2])
ax.axvline(0, linestyle='--', color='k')
ax.axhline(0, linestyle=':', color='k')
ax.legend(['no saccade','with saccade'])
ax.set_ylabel('L-R WSA [°]')
plt.tight_layout()
ax.xaxis.set_visible(False)
ax.set_yticks([-30,0,30,60])
ax.set_xlim(PLOT_TIMEWINDOW)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['right'].set_visible(False)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, hidey=False, sizey=0.0, barwidth=3)
<__main__.AnchoredScaleBar at 0x23588923988>
preSacc = 0.8 #s
postSacc = 0.5 #s
h=0.1 #m
v=1.5 #m/s
t0 = -h/v*1/np.tan(76.5/2) #76.5 is maximal angle (=34[px]*2.25[°/px])
theta = lambda t: 2*np.arctan(-h/v/(t+t0))*180/np.pi #in °
t_loom = np.arange(0.0, -preSacc, -1/50.)[::-1]
t_after = np.arange(0, postSacc, 1/50)
t_all = np.append(t_loom,t_after)
f,ax=plt.subplots(1,1)
ax.plot(t_all, np.append(theta(t_loom),np.full_like(t_after, 2.25)), color=[0.5]*3)
ax.axvline(0, linestyle='--', color='k')
plt.tight_layout()
ax.xaxis.set_visible(False)
ax.set_yticks([0,60])
ax.set_xlim(PLOT_TIMEWINDOW)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['right'].set_visible(False)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, hidey=False, sizey=0.0, barwidth=3)
plt.tight_layout()
# LOAD ALL SACCADE TRIALS OF ALL FLIES
anlyz=[None]*2
anlyz[0] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_saccading.neomat'))
anlyz[1] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_NoPol__CWTdetection_saccading.neomat'))
# anlyz[2] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight',
# 'saccLeft_DePol__CWTdetection_saccading.neomat'))
saSacc = sum(anlyz)
saSacc.average_flywise()
saSacc.flyID
{0: 'fly1_WTCS_20180725_cell1',
1: 'fly4_X21_20200420_cell1',
2: 'fly1_U29xG06_20190213_cell1',
3: 'fly1_WTCS_20190726_cell1',
4: 'fly3_U29xG16_20200129_cell1',
5: 'fly1_X20a_20200316_cell1',
6: 'fly1_G06xU29_20191104_cell1',
7: 'fly1_U29xG16_20191209_cell2',
8: 'fly3_X21_20200415_cell1',
9: 'fly3_x21_20200417_cell1',
10: 'fly2_X20_20210217_cell1'}
# LOAD ALL STRAIGHT TRIALS OF ALL FLIES
# anlyz=[None,None, None]
anlyz=[None,None]
anlyz[0] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_straight.neomat'))
anlyz[1] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_NoPol__CWTdetection_straight.neomat'))
# anlyz[2] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight',
# 'saccLeft_DePol__CWTdetection_straight.neomat'))
saStrg = sum(anlyz)
saStrg.average_flywise()
saStrg.flyID
{0: 'fly1_WTCS_20180725_cell1',
1: 'fly4_X21_20200420_cell1',
2: 'fly1_U29xG06_20190213_cell1',
3: 'fly1_WTCS_20190726_cell1',
4: 'fly3_U29xG16_20200129_cell1',
5: 'fly1_X20a_20200316_cell1',
6: 'fly1_G06xU29_20191104_cell1',
7: 'fly1_U29xG16_20191209_cell2',
8: 'fly3_X21_20200415_cell1',
9: 'fly3_x21_20200417_cell1',
10: 'fly2_X20_20210217_cell1'}
L-R WBA of saccades and non-saccade trials
preSacc = 0.7 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc.samplingRate), int(maxLoom+postSacc*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate-preSacc#-saSacc[i].preStim
f,ax=plt.subplots(1,1)
d = saStrg.get_avg_data('wingDiff')
s = saStrg.get_std_data('wingDiff')
# convert from radians to degrees
d = {k:v*180./np.pi for k,v in d.items()}
s = {k:v*180./np.pi for k,v in s.items()}
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]], color=plotColors[0])
ax.fill_between(timeSlice, (d[5]+s[5])[timeLims[0]:timeLims[1]], (d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc.get_avg_data('wingDiff')
s = saSacc.get_std_data('wingDiff')
# convert from radians to degrees
d = {k:v*180./np.pi for k,v in d.items()}
s = {k:v*180./np.pi for k,v in s.items()}
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]], color=plotColors[2])
ax.fill_between(timeSlice, (d[5]+s[5])[timeLims[0]:timeLims[1]], (d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
ax.set_ylabel('L-R WBA [°]')
ax.set_xlabel('time [s]')
ax.legend(['unresponsive','saccades'])
ax.axvline(0., color='k', linestyle=':') #marks ends of looming
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=5., labely='5 °', barwidth=3)
plt.tight_layout()
Run section before figure 2A first to load all the data
preSacc = 0.7 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc.samplingRate), int(maxLoom+postSacc*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate-preSacc#-saSacc[i].preStim
f,ax=plt.subplots(1,1)
d = saStrg.get_avg_data('ephys')
s = saStrg.get_std_data('ephys')
# conversion from raw amplifier output to mV
d = {k:v*100. for k,v in d.items()}
s = {k:v*100. for k,v in s.items()}
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]], color=plotColors[0])
ax.fill_between(timeSlice, (d[5]+s[5])[timeLims[0]:timeLims[1]], (d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc.get_avg_data('ephys')
s = saSacc.get_std_data('ephys')
# conversion from raw amplifier output to mV
d = {k:v*100. for k,v in d.items()}
s = {k:v*100. for k,v in s.items()}
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]], color=plotColors[2])
ax.fill_between(timeSlice, (d[5]+s[5])[timeLims[0]:timeLims[1]], (d[5]-s[5])[timeLims[0]:timeLims[1]], color=plotColors[2], alpha=plotAlphas[2])
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.legend(['straight','saccades'])
ax.axvline(0., color='k', linestyle=':') # marks end of looming
ax.set_xlim(PLOT_TIMEWINDOW)
ax.set_ylim([-3.5,4.5])
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=1., labely='1 mV', barwidth=3)
plt.tight_layout()
loads data and calculates amplitudes of membrane polarizations and L-R WSA for spontaneous saccades tp the left for setting line colors
folder = os.path.join(DATAPATH, 'spontaneousSaccades','CWTdetection')
#dirlist = os.lsdir(folder)
flist = glob.glob(os.path.join(folder,'*.neomat'))
remove = [i for i,f in enumerate(flist) if '20191025' in f] #too noisy, not used in other plots either
assert(len(remove)==1)
del flist[remove[0]]
flist
['D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_g06xu29_20191104_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_u29xg06_20190213_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_u29xg16_20191209_cell2_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_u29xg16_20200129_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_wtcs_20180725_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_wtcs_20190726_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_wtcs_20191106_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly1_x20a_20200316_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly2_x20_20210217_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly3_u29xg16_20200129_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly3_x21_20200415_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly3_x21_20200417_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly4_x21_20200420_cell1_looming_wb_spontSaccades_cwtDetection.neomat']
allSaccLeftMeanEphys = {}
allSaccLeftMeanWBA = {}
allSaccRightMeanEphys = {}
allSaccRightMeanWBA = {}
for fn in flist:
# load file into appropriate data structure (minimal subclass of BaseClass that keeps inter-trial intervals)
da = DummyAnalyzer.retrieveFromNeomatfile(fn)
#in some instances the ephys data has been saved as *100 (like used for plotting), so put all on same scale
if max(abs(da.saccLeftMeanEphys)) >0.1:
da.saccLeftMeanEphys /= 100.
da.saccRightMeanEphys /= 100.
# apply resampling where necessary
if abs(da.samplingRate -10000) >10: #not >0 because original double from MatLab was read with sligh deviation by Python (>1 should be sufficient)
#print('resampling ',fn)
da.saccLeftMeanEphys = signal.resample(da.saccLeftMeanEphys, int(da.saccLeftMeanEphys.size/da.samplingRate*10000.))
da.saccRightMeanEphys = signal.resample(da.saccRightMeanEphys, int(da.saccRightMeanEphys.size/da.samplingRate*10000.))
da.saccLeftMeanWBA = signal.resample(da.saccLeftMeanWBA, int(da.saccLeftMeanWBA.size/da.samplingRate*10000.))
da.saccRightMeanWBA = signal.resample(da.saccRightMeanWBA, int(da.saccRightMeanWBA.size/da.samplingRate*10000.))
# keep results in dict
allSaccLeftMeanEphys[da.flyID[0]] = da.saccLeftMeanEphys
allSaccRightMeanEphys[da.flyID[0]] = da.saccRightMeanEphys
# convert L-R WSA (=WBA) to degrees
allSaccLeftMeanWBA[da.flyID[0]] = da.saccLeftMeanWBA*180./np.pi
allSaccRightMeanWBA[da.flyID[0]] = da.saccRightMeanWBA*180./np.pi
plotlist = list(allSaccLeftMeanEphys.keys())[:]
#determine amplitudes of membrane polarizations and L-R WSA to attribute colors to line plots
polAmp=[] #polarizationAmplitude
polAmp_dict={}
saccAmp=[] #saccadeAmplitude
saccAmp_dict={}
for flyid in plotlist:
d=allSaccLeftMeanEphys[flyid]
d = ndimage.gaussian_filter1d(d, 60)
ix = int(np.argmax(signal.medfilt(np.abs(d[int(0.15*10000):]), 101))+0.15*10000)
polAmp.append(d[ix])
polAmp_dict[flyid] = d[ix]
d=allSaccLeftMeanWBA[flyid]
d = ndimage.gaussian_filter1d(d, 60)
ix = int(np.argmax(signal.medfilt(np.abs(d[int(0.15*10000):]), 101))+0.15*10000)
saccAmp.append(d[ix])
saccAmp_dict[flyid] = d[ix]
polAmp = np.array(polAmp)
saccAmp = np.array(saccAmp)
Plot L-R WBA of spontaneous saccades to the left, meaning R>L and L-R<0
assert(all([list(polAmp_dict.keys())[i]==plotlist[i] for i in range(len(plotlist))])) #make sure order is correct
polAmp = np.array([val for _,val in polAmp_dict.items()])
plotcolors = plt.cm.viridis((polAmp-min(polAmp))/(max(polAmp)-min(polAmp)))
preOnset = 0.15#s , fixed during saccade detection routine (bould could be set differently any time by reassigning saccLeftMeanWBA etc)
samplingRate = 10000. # or set to da.samplingRate
f,ax=plt.subplots(1,1)
for i,k in enumerate(plotlist):
ax.plot(np.arange(allSaccLeftMeanWBA[k].size)/samplingRate-preOnset, allSaccLeftMeanWBA[k],
color=plotcolors[i], linewidth=0.8)
ax.set_title('sacc. to left WSA')
ax.set_ylabel('$\Delta$ WSA [°]')
ax.set_xlabel('time from onset [s]')
ax.axvline(0, linestyle=':', color='k')
ax.set_xlim([-preOnset,0.51])
# add a line for average
# pad to longest to make dimensions match
padLen = max([v.size for k,v in allSaccLeftMeanWBA.items()])
allTraces = np.array([np.pad(v,(0,padLen-v.size), 'constant', constant_values=np.nan) for _,v in allSaccLeftMeanWBA.items()])
ax.plot(np.arange(padLen)/10000.-0.15, np.nanmean(allTraces,axis=0), 'k-', linewidth=2.5)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=5., labely='5 °', barwidth=3)
plt.tight_layout()
membrane potential corresponding to traces in Fig 2C
assert(all([list(polAmp_dict.keys())[i]==plotlist[i] for i in range(len(plotlist))])) #make sure order is correct
polAmp = np.array([val for _,val in polAmp_dict.items()])
plotcolors = plt.cm.viridis((polAmp-min(polAmp))/(max(polAmp)-min(polAmp)))
preOnset = 0.15#s , fixed during saccade detection routine (bould could be set differently any time by reassigning saccLeftMeanEphys etc)
samplingRate = 10000. # or set to da.samplingRate
f,ax=plt.subplots(1,1)
for i,k in enumerate(plotlist):
ax.plot(np.arange(allSaccLeftMeanEphys[k].size)/samplingRate-preOnset, allSaccLeftMeanEphys[k]*100.,
color=plotcolors[i], linewidth=1.5)
ax.set_title('sacc. to left membrane potential')
ax.set_ylabel('membrane potential [mV]')
ax.set_xlabel('time from onset [s]')
ax.axvline(0, linestyle=':', color='k')
ax.set_xlim([-preOnset,0.51])
# add a line for average
# pad to longest to make dimensions match
padLen = max([v.size for k,v in allSaccLeftMeanEphys.items()])
allTraces = np.array([np.pad(v,(0,padLen-v.size), 'constant', constant_values=np.nan) for _,v in allSaccLeftMeanEphys.items()])
ax.plot(np.arange(padLen)/10000.-0.15, np.nanmean(allTraces,axis=0)*100., 'k-', linewidth=2.5)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2., labely='2 mV', barwidth=3)
plt.tight_layout()
saSacc = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_saccading.neomat'))
saSacc.average_flywise()
saStrg = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_straight.neomat'))
saStrg.average_flywise()
# prettier, clean plotting, membrane potential
# assert that we have equal datasets
strgSet = set([v.lower() for _,v in saStrg.flyID.items()])
saccSet = set([v.lower() for _,v in saSacc.flyID.items()])
assert(strgSet == saccSet)
f,ax=plt.subplots(1,1)
preSacc = 0.7 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc.samplingRate), int(maxLoom+postSacc*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate-preSacc#-saSacc[i].preStim
d = saStrg.get_avg_data('ephys')
s = saStrg.get_std_data('ephys')
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]]*100.)
ax.fill_between(timeSlice, 100.*(d[5]+s[5])[timeLims[0]:timeLims[1]], 100.*(d[5]-s[5])[timeLims[0]:timeLims[1]], alpha=0.3)
d = saSacc.get_avg_data('ephys')
s = saSacc.get_std_data('ephys')
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]]*100.)
ax.fill_between(timeSlice, 100.*(d[5]+s[5])[timeLims[0]:timeLims[1]], 100.*(d[5]-s[5])[timeLims[0]:timeLims[1]], alpha=0.3)
ax.set_title('saccades to left')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preSacc, postSacc])
ax.axvline(0, linestyle=':', color='k')
ax.legend(['unresponsive','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
ax.set_ylim([-3.5,4.5])
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=1., labely='1 mV', barwidth=3)
plt.tight_layout()
saSacc = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_NoPol__CWTdetection_saccading.neomat'))
saSacc.average_flywise()
saStrg = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH, 'Collections_saccade_straight',
'saccLeft_NoPol__CWTdetection_straight.neomat'))
saStrg.average_flywise()
# membrane potential, multiply by factor of 100 to convert amplifier output to mV
# assert that we have equal datasets
strgSet = set([v for _,v in saStrg.flyID.items()])
saccSet = set([v for _,v in saSacc.flyID.items()])
assert(strgSet == saccSet)
f,ax=plt.subplots(1,1)
preSacc = 0.7 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc.samplingRate), int(maxLoom+postSacc*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate-preSacc
d = saStrg.get_avg_data('ephys')
s = saStrg.get_std_data('ephys')
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[5]+s[5])[timeLims[0]:timeLims[1]], 100.*(d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc.get_avg_data('ephys')
s = saSacc.get_std_data('ephys')
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[5]+s[5])[timeLims[0]:timeLims[1]], 100.*(d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
# ax.set_title('saccades to left')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle=':', color='k')
ax.legend(['unresponsive','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
ax.set_ylim([-3.5,4.5])
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=1., labely='1 mV', barwidth=3)
plt.tight_layout()
Plotting membrane potential during spontaneous and escape saccades side-by-side for the hyperpolarizing group
anlyz = []
flist= os.listdir(os.path.join(DATAPATH,'spontaneousSaccades','CWTdetection'))
#this simplifies in other folder structures
for fly in saccLeft_HyPol:
#print(fly)
fn = fly[0]
if 'vonBettina' in fn:
date = fn.split('/')[1][:8]
date = date[-4:]+date[2:4]+date[:2]
else:
date = ''.join(fn.split('_'))
name = '_'.join(fly[1].lower().split('_')[:2])
name = name+'_'+date
print(name)
fl = [f for f in flist if name in f]
fl = [f for f in fl if f.endswith('.neomat')]
assert(len(fl)==1)
a = DummyAnalyzer.retrieveFromNeomatfile(os.path.join(DATAPATH, 'spontaneousSaccades','CWTdetection', fl[0]))
a.group_indices[a.group_indices!=0]=0
helperLibrary.resampleAnalyzer(a) #wills kip if sampling rate is already 10kHz
anlyz.append(a)
fly1_wtcs_20180725 fly1_u29xg06_20190213 fly1_wtcs_20190726 fly3_u29xg16_20200129 fly1_x20a_20200316 fly4_x21_20200420
def extract_saccade_LR_ephys_WBA(anlyz, **kwargs):
allLeftEphys = []
allLeftWBA = []
allRightEphys = []
allRightWBA = []
preOnset_t = kwargs.get('preOnset',0.15)
postOnset_t = kwargs.get('afterOnset',0.5)
for el in range(len(anlyz)):
if not isinstance(anlyz.saccades[el]['peakTime'], (list,np.ndarray)):
anlyz.saccades[el] = {k:np.array([v]) for k,v in anlyz.saccades[el].items()}
if len(anlyz.saccades[el]['peakTime'])==0: continue
v = anlyz.saccades[el]
actualOnset= (v['CWTpeaks']*0.02*anlyz.samplingRate).astype(int)
saccAmp = v['saccAmp']
rightSaccEphys = []
leftSaccEphys=[]
leftSaccWBA=[]
rightSaccWBA=[]
chunk = int(postOnset_t*anlyz.samplingRate)
preOnset = int(preOnset_t*anlyz.samplingRate)
ephys = anlyz.segments[el].data[0]*100.#ephys data
wba = anlyz.segments[el].data[-1]
# CONVERT
wba = wba*180./np.pi
for j in range(saccAmp.size):
o=actualOnset[j]
a=saccAmp[j]
p= (v['peakTime']*0.02*anlyz.samplingRate).astype(int)[j]
if a<0:
leftSaccEphys.append(ephys[o-preOnset:o+chunk])
leftSaccWBA.append(wba[o-preOnset:o+chunk])
d = wba[o-preOnset:o+chunk]
xval = np.arange(-preOnset, d.size-preOnset)/anlyz.samplingRate
d = ephys[o-preOnset:o+chunk]
else:
rightSaccEphys.append(ephys[o-preOnset:o+chunk])
rightSaccWBA.append(wba[o-preOnset:o+chunk])
d = wba[o-preOnset:o+chunk]
xval = np.arange(-preOnset,d.size-preOnset)/anlyz.samplingRate
d = ephys[o-preOnset:o+chunk]
allLeftEphys.extend(leftSaccEphys)
allLeftWBA.extend(leftSaccWBA)
allRightEphys.extend(rightSaccEphys)
allRightWBA.extend(rightSaccWBA)
return allLeftEphys, allLeftWBA, allRightEphys, allRightWBA
# AVERAGE OVER SACCADES AND STD
allLeftEphys = []
allLeftWBA = []
preOnset = 0.15#s
for a in anlyz:
leftEphys,leftWBA,_,_ = extract_saccade_LR_ephys_WBA(a, preOnset=preOnset)
allLeftEphys.extend(leftEphys)
allLeftWBA.extend(leftWBA)
f,ax=plt.subplots(1,1)
maxLen = int(np.max([len(y) for y in allLeftEphys]))
xval = np.arange(maxLen)/anlyz[0].samplingRate-preOnset
m = np.nanmean(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan) for y in allLeftEphys]), axis=0)
s = np.nanstd(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan) for y in allLeftEphys]), axis=0)
# s /= np.sqrt(len(anlyz))
ax.plot(xval, m, color='g')
ax.fill_between(xval, m+s,m-s, color='g', alpha=0.4)
# ax.set_title('spontaneous saccades to left')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preOnset, 0.5])
ax.axvline(0, linestyle=':', color='k')
# plt.tight_layout()
<matplotlib.lines.Line2D at 0x235a00d3b08>
# AVERAGE SACCADE, average flywise and proper s.e.m
# ephys
allLeftEphys = []
allLeftEphys_std = []
preOnset = 0.15#s
for a in anlyz:
leftEphys,leftWBA,_,_ = extract_saccade_LR_ephys_WBA(a, preOnset=preOnset)
maxLen = max([len(y) for y in leftEphys])
meanLeftEphys = np.nanmean(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan)
for y in leftEphys]), axis=0)
stdLeftEphys = np.nanstd(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan)
for y in leftEphys]), axis=0)
allLeftEphys.append(meanLeftEphys)
allLeftEphys_std.append(stdLeftEphys)
f,ax=plt.subplots(1,1)
maxLen = int(np.max([len(y) for y in allLeftEphys]))
xval = np.arange(maxLen)/anlyz[0].samplingRate-preOnset
m = np.nanmean(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan) for y in allLeftEphys]), axis=0)
s = np.sqrt(np.sum(np.array([np.pad(y,(0,maxLen-len(y)), mode='constant', constant_values=np.nan)**2.
for y in allLeftEphys_std]),axis=0) / len(anlyz)**2)
ax.plot(xval, m, color='#212178')
ax.fill_between(xval, m+s,m-s, color='#212178', alpha=0.4)
ax.set_title('spontaneous saccades to left')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preOnset, 0.5])
ax.axvline(0, linestyle=':', color='k')
<matplotlib.lines.Line2D at 0x23598191b48>
saSacc = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_saccading.neomat'))
saSacc.flyID
{0: 'fly1_WTCS_20180725_cell1',
1: 'fly4_X21_20200420_cell1',
2: 'fly1_U29xG06_20190213_cell1',
3: 'fly1_WTCS_20190726_cell1',
4: 'fly3_U29xG16_20200129_cell1',
5: 'fly1_X20a_20200316_cell1'}
# saSacc.average_by_stimSignal()
saSacc.average_flywise()
preSacc = 0.2 #s
postSacc = 0.5 #s
maxLoom = (saSacc.preStim+saSacc.stimulusPhases_inframes[1]*saSacc.meanFramePeriod)*saSacc.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc.samplingRate), int(maxLoom+postSacc*saSacc.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc.samplingRate-preSacc#-saSacc[i].preStim
d = saSacc.get_avg_data('ephys')
s = saSacc.get_std_data('ephys')
ax.plot(timeSlice, d[5][timeLims[0]:timeLims[1]]*100., color=plotColors[2], linewidth=1.5)
ax.fill_between(timeSlice, 100.*(d[5]+s[5])[timeLims[0]:timeLims[1]], 100.*(d[5]-s[5])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
<matplotlib.collections.PolyCollection at 0x235881f1a48>
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=1., labely='1 mV', barwidth=3)
<__main__.AnchoredScaleBar at 0x229dc235b48>
# load new subsets
saSacc_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_straight.neomat'))
saSacc_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_saccading.neomat'))
saStrg_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_straight.neomat'))
helperLibrary.resampleAnalyzer(saSacc_Rhyp)
helperLibrary.resampleAnalyzer(saStrg_Rhyp)
helperLibrary.resampleAnalyzer(saSacc_Rnop)
helperLibrary.resampleAnalyzer(saStrg_Rnop)
saStrg_Rhyp.clean_flyIDs()
saStrg_Rhyp.reallocate_fly_indices()
saStrg_Rnop.clean_flyIDs()
saStrg_Rnop.reallocate_fly_indices()
saSacc_Rhyp.clean_flyIDs()
saSacc_Rhyp.reallocate_fly_indices()
saSacc_Rnop.clean_flyIDs()
saSacc_Rnop.reallocate_fly_indices()
saSacc_all = copy.copy(saSacc_Rhyp)
saStrg_all = copy.copy(saStrg_Rhyp)
#saFly = copy.copy(saFly_Rhyp)
#saRest = copy.copy(saRest_Rhyp)
saSacc_all = saSacc_all + saSacc_Rnop
saStrg_all = saStrg_all + saStrg_Rnop
#saFly = saFly + saFly_Rnop
#saRest = saRest + saRest_Rnop
saSacc_all.baseline_correct_data(baseline_startpoint=int((saSacc_all.preStim+0.15)*saSacc_all.samplingRate))
saStrg_all.baseline_correct_data(baseline_startpoint=int((saStrg_all.preStim+0.15)*saStrg_all.samplingRate))
saSacc_all.average_by_stimSignal()
saStrg_all.average_by_stimSignal()
# prettier, clean plotting, membrane potential
# assert that we have equal datasets
strgSet = set([v.lower() for _,v in saStrg_all.flyID.items()])
# restSet = set([v.lower() for _,v in saRest.flyID.items()])
saccSet = set([v.lower() for _,v in saSacc_all.flyID.items()])
''' MASSIVE NO-GO FOR NOW; BECAUSE I NEED TO FIX THE FLYID THING'''
#assert(strgSet == restSet)
#assert(saccSet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.6 #s
postSacc = 0.5 #s
maxLoom = (saSacc_all.preStim+saSacc_all.stimulusPhases_inframes[1]*saSacc_all.meanFramePeriod)*saSacc_all.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_all.samplingRate), int(maxLoom+postSacc*saSacc_all.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_all.samplingRate-preSacc#-saSacc_all[i].preStim
# d = saRest.get_avg_data('wingDiff')
# s = saRest.get_std_data('wingDiff')
# ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*180./np.pi, color=plotColors[0])
# ax.fill_between(timeSlice, 180./np.pi*(d[2]+s[2])[timeLims[0]:timeLims[1]], 180./np.pi*(d[2]-s[2])[timeLims[0]:timeLims[1]],
# color=plotColors[0], alpha=plotAlphas[0])
d = saStrg_all.get_avg_data('wingDiff')
s = saStrg_all.get_std_data('wingDiff')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*180./np.pi, color=plotColors[0])
ax.fill_between(timeSlice, 180./np.pi*(d[2]+s[2])[timeLims[0]:timeLims[1]], 180./np.pi*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc_all.get_avg_data('wingDiff')
s = saSacc_all.get_std_data('wingDiff')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*180./np.pi, color=plotColors[2])
ax.fill_between(timeSlice, 180./np.pi*(d[2]+s[2])[timeLims[0]:timeLims[1]], 180./np.pi*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preSacc, postSacc])
ax.axvline(0, linestyle=':', color='k')
# ax.legend(['rest', 'no saccades','saccades'])
ax.legend(['no saccades','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
#f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=10.0, labely='10 °', barwidth=3)
plt.tight_layout()
# membrane potential
# assert that we have equal datasets
strgSet = set([v.lower() for _,v in saStrg_all.flyID.items()])
saccSet = set([v.lower() for _,v in saSacc_all.flyID.items()])
''' MASSIVE NO-GO FOR NOW; BECAUSE I NEED TO FIX THE FLYID THING'''
#assert(strgSet == restSet)
#assert(saccSet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.6 #s
postSacc = 0.5 #s
maxLoom = (saSacc_all.preStim+saSacc_all.stimulusPhases_inframes[1]*saSacc_all.meanFramePeriod)*saSacc_all.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_all.samplingRate), int(maxLoom+postSacc*saSacc_all.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_all.samplingRate-preSacc#-saSacc_all[i].preStim
d = saStrg_all.get_avg_data('ephys')
s = saStrg_all.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[0])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc_all.get_avg_data('ephys')
s = saSacc_all.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preSacc, postSacc])
ax.axvline(0, linestyle=':', color='k')
ax.legend(['no saccades','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
ax.set_ylim([-4.8,6.5]) #avgSEM
# ax.set_ylim([-6.,8.35]) #avgSTD-byStimSig
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
The left part of this panel is a purely artistic schematic drawing. The right part is data-derived, though.
allflies = saccLeft_HyPol[:]
allflies.extend(saccLeft_DePol)
allflies.extend(saccLeft_NoPol)
# LOAD
saFly = []
for root,fly in allflies:
a = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH, root, f'{fly}looming_wb__flying.neomat'))
helperLibrary.resampleAnalyzer(a)# only does something where necessary
a.average_flywise()
saFly.append(a)
# other way, using collections (originally established only later), previous cell kept for demonstration
anlyz = [None, None]
anlyz[0] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__flying.neomat'))
anlyz[1] = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__flying.neomat'))
saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile()
saSacc_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_saccading.neomat'))
# LOAD
saRest = []
for root,fly in allflies:
a = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH, root, f'{fly}looming_wb__resting.neomat'))
helperLibrary.resampleAnalyzer(a)# only does something where necessary
a.average_flywise()
saRest.append(a)
# extract a mean trace overa ll flies (without having to combine all datasets, which is tedious)
# -> not necessary, only for confirming other results
stimID = 2
loomLeft=[]
legEntry=[]
for a in saFly:
d = a.get_avg_data('ephys')[stimID]*100.# conversion factor to mV
d = ndimage.gaussian_filter1d(d,5)
loomLeft.append(d)
legEntry.append(a.flyID[0])
minLen = min([len(d) for d in loomLeft]) # cut to shortest instance (which still covers all the stimulus anyhow)
loomLeft = np.array([d[:minLen] for d in loomLeft])
time = a.timeax
meanFlight = np.mean(loomLeft,axis=0)
before = []
after = []
deltas = []
beforerange = np.arange(7000,8000) #stretch of data (indices) beofre maximum looming
afterrange = np.arange(11000,12000) #stretch of data (indices) after maximum looming
f,ax=plt.subplots(1,1, figsize=(3,5))
rng = np.random.default_rng() # for creating slight x-offset so points don't cover one another
for a in saFly:
fid = a.flyID[0]
print(fid)
d = a.get_avg_data('ephys')[stimID]*100.
d = ndimage.gaussian_filter1d(d,5)
b = np.mean(d[beforerange])
before.append(b)
a = np.mean(d[afterrange])
after.append(a)
deltas.append(a-b)
ax.plot(0+rng.uniform(low=-0.001,high=0.001), a-b,'o', markersize=10)
meanDelta = np.mean(deltas)
#stdDelta = np.std(deltas)
ax.plot(meanDelta, 'ks', markersize=11)
#ax.plot([0,0], [meanDelta-stdDelta, meanDelta+stdDelta], 'k-')
#prettify
ax.set_xlim([-0.01,0.01])
ax.set_ylabel('$\mu_{after}-\mu_{before}$', fontsize=20)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.xaxis.set_ticks([])
ax.yaxis.set_tick_params(labelsize=15)
ax.yaxis.set_ticks([-4,-2,0,2])
plt.tight_layout()
fly1_WTCS_20180725_cell1 fly1_WTCS_20190213_cell1 fly1_WTCS_20190726_cell1 fly3_U29xG16_20200129_cell1 fly1_X20a_20200316_cell1 fly4_X21_20200420_cell1 fly1_WTCS_20191106_cell1 fly1_U29xG16_20200129_cell1 fly1_G06xU29_20191104_cell1 fly1_U29xG16_20191209_cell2 fly3_X21_20200415_cell1 fly3_x21_20200417_cell1 fly2_X20_20210217_cell1
# CHANGE COLORS AFTER KNOWING THE MEAN
# alternative way to retrieve diffs from plot
#diffs = np.array([ch.get_ydata()[0] for ch in ax[1].get_children() if isinstance(ch,plt.Line2D)]).flatten()
dots = [ch for ch in ax.get_children() if isinstance(ch, plt.Line2D)]
idxAboveMean = np.where(deltas>meanDelta)[0]
dotColors = plt.cm.YlGn(np.linspace(0.25,0.75,len(idxAboveMean)))
# dotColors = plt.cm.Blues(np.linspace(0.3,0.6,len(idxAboveMean)))
for ix,col in zip(idxAboveMean, dotColors):
dots[ix].set_color(col)
idxBelowMean = np.where(deltas<=meanDelta)[0]
dotColors = plt.cm.Purples(np.linspace(0.45,0.7,len(idxBelowMean)))
# dotColors = plt.cm.Oranges(np.linspace(0.3,0.65,len(idxBelowMean)))
for ix,col in zip(idxBelowMean, dotColors):
dots[ix].set_color(col)
Load datasets although we did that before in section 5.1, but it shoudln't be necessary to run that first.
saSacc_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_straight.neomat'))
saSacc_Rhyp.clean_flyIDs()
saStrg_Rhyp.clean_flyIDs()
saSacc_Rhyp.reallocate_fly_indices()
saStrg_Rhyp.reallocate_fly_indices()
saSacc_Rhyp.baseline_correct_data(baseline_startpoint=int((saSacc_Rhyp.preStim+0.15)*saSacc_Rhyp.samplingRate))
saStrg_Rhyp.baseline_correct_data(baseline_startpoint=int((saStrg_Rhyp.preStim+0.15)*saStrg_Rhyp.samplingRate))
saSacc_Rhyp.average_by_stimSignal()
saStrg_Rhyp.average_by_stimSignal()
# membrane potential
# assert that we have equal datasets
strgSet = set([v for _,v in saStrg_Rhyp.flyID.items()])
saccSet = set([v for _,v in saSacc_Rhyp.flyID.items()])
''' MASSIVE NO-GO FOR NOW; BECAUSE I NEED TO FIX THE FLYID THING'''
#assert(strgSet == restSet)
#assert(saccSet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.6 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rhyp.preStim+saSacc_Rhyp.stimulusPhases_inframes[1]*saSacc_Rhyp.meanFramePeriod)*saSacc_Rhyp.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rhyp.samplingRate), int(maxLoom+postSacc*saSacc_Rhyp.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rhyp.samplingRate-preSacc#-saSacc_Rhyp[i].preStim
d = saStrg_Rhyp.get_avg_data('ephys')
s = saStrg_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[0])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc_Rhyp.get_avg_data('ephys')
s = saSacc_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preSacc, postSacc])
ax.axvline(0, linestyle=':', color='k')
ax.legend(['no saccades','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
# ax.set_ylim([-4.8,6.5]) #avgSEM
ax.set_ylim([-6.,8.35]) #avgSTD-byStimSig
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
Load datasets although we did that before in section 5.1, but it shoudln't be necessary to run that first.
saSacc_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_saccading.neomat'))
saStrg_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_straight.neomat'))
saSacc_Rnop.clean_flyIDs()
saStrg_Rnop.clean_flyIDs()
saSacc_Rnop.reallocate_fly_indices()
saStrg_Rnop.reallocate_fly_indices()
saSacc_Rnop.baseline_correct_data(baseline_startpoint=int((saSacc_Rnop.preStim+0.15)*saSacc_Rnop.samplingRate))
saStrg_Rnop.baseline_correct_data(baseline_startpoint=int((saStrg_Rnop.preStim+0.15)*saStrg_Rnop.samplingRate))
saSacc_Rnop.average_by_stimSignal()
saStrg_Rnop.average_by_stimSignal()
# membrane potential
# assert that we have equal datasets
strgSet = set([v for _,v in saStrg_Rnop.flyID.items()])
saccSet = set([v for _,v in saSacc_Rnop.flyID.items()])
''' MASSIVE NO-GO FOR NOW; BECAUSE I NEED TO FIX THE FLYID THING'''
#assert(strgSet == restSet)
#assert(saccSet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.6 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rnop.preStim+saSacc_Rnop.stimulusPhases_inframes[1]*saSacc_Rnop.meanFramePeriod)*saSacc_Rnop.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rnop.samplingRate), int(maxLoom+postSacc*saSacc_Rnop.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rnop.samplingRate-preSacc#-saSacc_Rnop[i].preStim
d = saStrg_Rnop.get_avg_data('ephys')
s = saStrg_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[0])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[0], alpha=plotAlphas[0])
d = saSacc_Rnop.get_avg_data('ephys')
s = saSacc_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[2][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[2]+s[2])[timeLims[0]:timeLims[1]], 100.*(d[2]-s[2])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim([-preSacc, postSacc])
ax.axvline(0, linestyle=':', color='k')
ax.legend(['no saccades','saccades'])
ax.set_xlim(PLOT_TIMEWINDOW)
# ax.set_ylim([-4.8,6.5]) #avgSEM
ax.set_ylim([-6.,8.35]) #avgSTD-byStimSig
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
# LOAD DATASETS JUST TO GET THE LIST OF FLIES ASSOCIATED WITH EACH GROUP
# IF you have an alternative way, go ahead. Like using the ids straight from Fig 3C or D&E instead
saFly_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__flying.neomat'))
saFly_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__flying.neomat'))
saFly_Rhyp.clean_flyIDs()
saFly_Rnop.clean_flyIDs()
saFly_Rhyp.reallocate_fly_indices()
saFly_Rnop.reallocate_fly_indices()
hyPolSet = set([flyid.lower() for _,flyid in saFly_Rhyp.flyID.items()])
noPolSet = set([flyid.lower() for _,flyid in saFly_Rnop.flyID.items()])
assert(noPolSet.isdisjoint(hyPolSet))
# in this figure _variableX_noPol refers to those from the NP_contra group
# and _variableX_hyPol refers to those from the HP_contra group
allSaccLeftMeanEphys_noPol = {}
allSaccLeftMeanWBA_noPol = {}
allSaccRightMeanEphys_noPol = {}
allSaccRightMeanWBA_noPol = {}
allSaccLeftMeanEphys_hyPol = {}
allSaccLeftMeanWBA_hyPol = {}
allSaccRightMeanEphys_hyPol = {}
allSaccRightMeanWBA_hyPol = {}
# load NP_contra spontaneous saccade data, relying on existing data structure
for flyID in noPolSet:
print(flyID)
fn = os.path.join(DATAPATH, 'spontaneousSaccades','CWTdetection',
f'{flyID}_looming_wb_spontSaccades_cwtDetection.neomat')
da = DummyAnalyzer.retrieveFromNeomatfile(fn)
#assert (da.flyID[0]==flyID.lower())
if max(abs(da.saccLeftMeanEphys)) >0.1:
da.saccLeftMeanEphys /= 100.
da.saccRightMeanEphys /= 100.
if abs(da.samplingRate -10000) >10:
da.saccLeftMeanEphys = signal.resample(da.saccLeftMeanEphys, int(da.saccLeftMeanEphys.size/da.samplingRate*10000.))
da.saccRightMeanEphys = signal.resample(da.saccRightMeanEphys, int(da.saccRightMeanEphys.size/da.samplingRate*10000.))
da.saccLeftMeanWBA = signal.resample(da.saccLeftMeanWBA, int(da.saccLeftMeanWBA.size/da.samplingRate*10000.))
da.saccRightMeanWBA = signal.resample(da.saccRightMeanWBA, int(da.saccRightMeanWBA.size/da.samplingRate*10000.))
allSaccLeftMeanEphys_noPol[da.flyID[0]] = da.saccLeftMeanEphys
allSaccRightMeanEphys_noPol[da.flyID[0]] = da.saccRightMeanEphys
allSaccLeftMeanWBA_noPol[da.flyID[0]] = da.saccLeftMeanWBA*180./np.pi
allSaccRightMeanWBA_noPol[da.flyID[0]] = da.saccRightMeanWBA*180./np.pi
# load HP_contra spontaneous saccade data, relying on existing data structure
for flyID in hyPolSet:
print(flyID)
fn = os.path.join(DATAPATH, 'spontaneousSaccades','CWTdetection',
f'{flyID}_looming_wb_spontSaccades_cwtDetection.neomat')
da = DummyAnalyzer.retrieveFromNeomatfile(fn)
#assert (da.flyID[0]==flyID.lower())
if max(abs(da.saccLeftMeanEphys)) >0.1:
da.saccLeftMeanEphys /= 100.
da.saccRightMeanEphys /= 100.
if abs(da.samplingRate -10000) >10:
da.saccLeftMeanEphys = signal.resample(da.saccLeftMeanEphys, int(da.saccLeftMeanEphys.size/da.samplingRate*10000.))
da.saccRightMeanEphys = signal.resample(da.saccRightMeanEphys, int(da.saccRightMeanEphys.size/da.samplingRate*10000.))
da.saccLeftMeanWBA = signal.resample(da.saccLeftMeanWBA, int(da.saccLeftMeanWBA.size/da.samplingRate*10000.))
da.saccRightMeanWBA = signal.resample(da.saccRightMeanWBA, int(da.saccRightMeanWBA.size/da.samplingRate*10000.))
allSaccLeftMeanEphys_hyPol[da.flyID[0]] = da.saccLeftMeanEphys
allSaccRightMeanEphys_hyPol[da.flyID[0]] = da.saccRightMeanEphys
allSaccLeftMeanWBA_hyPol[da.flyID[0]] = da.saccLeftMeanWBA*180./np.pi
allSaccRightMeanWBA_hyPol[da.flyID[0]] = da.saccRightMeanWBA*180./np.pi
fly1_u29xg16_20200129_cell1 fly3_x21_20200415_cell1 fly3_x21_20200417_cell1 fly1_g06xu29_20191104_cell1 fly1_u29xg16_20191209_cell2 fly1_x20a_20200316_cell1 fly3_u29xg16_20200129_cell1 fly1_wtcs_20191106_cell1 fly2_x20_20210217_cell1
Plot spontaneous saccades L-R WSA for subsets in one plot
f,ax = plt.subplots(1,1)
plotcolors = plt.cm.YlGn(np.linspace(0.3,0.6,len(allSaccRightMeanWBA_noPol)))
#plotcolors = plt.cm.YlOrBr(np.linspace(0.3,0.7,len(allSaccRightMeanWBA_noPol)))
i=0
for fid, data in allSaccRightMeanWBA_noPol.items():
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1.)
i+=1
# add average line
minLen = min([data.size for _,data in allSaccRightMeanWBA_noPol.items()])
alldata = np.array([data[:minLen] for _,data in allSaccRightMeanWBA_noPol.items()])
avg_Rnop = np.mean(alldata, axis=0)
del alldata, plotcolors
plotcolors = plt.cm.Purples(np.linspace(0.45,0.7,len(allSaccRightMeanWBA_hyPol)))
i=0
for fid, data in allSaccRightMeanWBA_hyPol.items():
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1.)
i+=1
# add average line
minLen = min([data.size for _,data in allSaccRightMeanWBA_hyPol.items()])
alldata = np.array([data[:minLen] for _,data in allSaccRightMeanWBA_hyPol.items()])
avg_Rhyp = np.mean(alldata, axis=0)
del alldata, plotcolors
# ax.set_title('sacc right WBA')
# ax.set_xlabel('time from onset [s]')
# compute and add global average (weighted average of the two group averages, number of members as weight)
minLen = min([avg_Rnop.size, avg_Rhyp.size])
avgAll = ((avg_Rnop[:minLen]*len(saFly_Rnop.flyID) + avg_Rhyp[:minLen]*len(saFly_Rhyp.flyID))
/ (len(saFly_Rnop.flyID) + len(saFly_Rhyp.flyID)) )
ax.plot(np.arange(avgAll.size)/10000-0.15, avgAll, color=[0.3,0.3,0.3], linewidth=2)
ax.plot(np.arange(avg_Rhyp.size)/10000.-0.15,avg_Rhyp, color=plt.cm.Purples(0.98), linewidth=2)
ax.plot(np.arange(avg_Rnop.size)/10000.-0.15,avg_Rnop, color=plt.cm.YlGn(0.85), linewidth=2.5)
ax.set_xlim([-0.15,0.48])
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=5., labely='5 °', barwidth=3)
plt.tight_layout()
# first plot NP_contra in green, then HP_contra in purple
# average lines per group and all flies total are plotted at the end so they end up on top
f,ax = plt.subplots(1,1)
plotcolors = plt.cm.YlGn(np.linspace(0.3,0.6,len(allSaccRightMeanEphys_noPol)))
#plotcolors = plt.cm.YlOrBr(np.linspace(0.3,0.7,len(allSaccRightMeanEphys_noPol)))
i=0
for fid, data in allSaccRightMeanEphys_noPol.items():
ax.plot(np.arange(data.size)/10000.-0.15, data*100., color=plotcolors[i], linewidth=1.)
i+=1
# add average line
minLen = min([data.size for _,data in allSaccRightMeanEphys_noPol.items()])
alldata = np.array([data[:minLen] for _,data in allSaccRightMeanEphys_noPol.items()])
avg_Rnop = np.mean(alldata, axis=0)
del alldata, plotcolors
plotcolors = plt.cm.Purples(np.linspace(0.45,0.7,len(allSaccRightMeanEphys_hyPol)))
i=0
for fid, data in allSaccRightMeanEphys_hyPol.items():
ax.plot(np.arange(data.size)/10000.-0.15, data*100., color=plotcolors[i], linewidth=1.)
i+=1
# add average line
minLen = min([data.size for _,data in allSaccRightMeanEphys_hyPol.items()])
alldata = np.array([data[:minLen] for _,data in allSaccRightMeanEphys_hyPol.items()])
avg_Rhyp = np.mean(alldata, axis=0)
del alldata, plotcolors
# compute and add global average (weighted average of the two group averages, number of members as weight)
minLen = min([avg_Rnop.size, avg_Rhyp.size])
avgAll = ((avg_Rnop[:minLen]*len(saFly_Rnop.flyID) + avg_Rhyp[:minLen]*len(saFly_Rhyp.flyID))
/ (len(saFly_Rnop.flyID) + len(saFly_Rhyp.flyID)) )
# plot average lines here in certain order for better visibility (they end up on top)
ax.plot(np.arange(avgAll.size)/10000-0.15, avgAll*100., color=[0.3,0.3,0.3], linewidth=3)
ax.plot(np.arange(avg_Rnop.size)/10000.-0.15,avg_Rnop*100., color=plt.cm.YlGn(0.85), linewidth=3.5)
ax.plot(np.arange(avg_Rhyp.size)/10000.-0.15,avg_Rhyp*100., color=plt.cm.Purples(0.98), linewidth=3)
ax.set_xlim([-0.15,0.48])
ax.set_title('sacc right Ephys')
ax.set_xlabel('time from onset [s]')
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=1.0, labely='1 mV', barwidth=3)
plt.tight_layout()
PLOT_Y_LIMS = np.array([-3, 8.])
PLOT_TIMEWINDOW = [-0.4,0.46]
PLOT_LEGEND = ['horizontal motion only','rest', 'saccades']
# simpler alternative
saFly_Lhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__flying.neomat'))
saRest_Lhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__resting.neomat'))
# saccading trials
saSacc_Lhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__CWTdetection_saccading.neomat'))
#saSacc_Lhyp.clean_flyIDs()
#saSacc_Lhyp.reallocate_fly_indices()
saSacc_Lhyp.flyID
{0: 'fly1_WTCS_20180725_cell1',
1: 'fly4_X21_20200420_cell1',
2: 'fly1_U29xG06_20190213_cell1',
3: 'fly1_WTCS_20190726_cell1',
4: 'fly3_U29xG16_20200129_cell1',
5: 'fly1_X20a_20200316_cell1'}
saRest_Lhyp.average_by_stimSignal()
saSacc_Lhyp.average_by_stimSignal()
saFly_Lhyp.average_by_stimSignal()
# prettier, clean plotting, membrane potential
stimID = 3
# assert that we have equal datasets
restSet = set([v.lower() for _,v in saRest_Lhyp.flyID.items()])
flySet = set([v.lower() for _,v in saSacc_Lhyp.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Lhyp.preStim+saSacc_Lhyp.stimulusPhases_inframes[1]*saSacc_Lhyp.meanFramePeriod)*saSacc_Lhyp.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Lhyp.samplingRate), int(maxLoom+postSacc*saSacc_Lhyp.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Lhyp.samplingRate-preSacc#-saSacc_Lhyp[i].preStim
d = saFly_Lhyp.get_avg_data('ephys')
s = saFly_Lhyp.get_std_data('ephys')
ax.plot(timeSlice, d[6][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[6]+s[6])[timeLims[0]:timeLims[1]], 100.*(d[6]-s[6])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Lhyp.get_avg_data('ephys')
s = saRest_Lhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Lhyp.get_avg_data('ephys')
s = saSacc_Lhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to left and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.set_xlim(PLOT_TIMEWINDOW)
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS)
#f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
# simpler alternative
saFly_Lhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__flying.neomat'))
saRest_Lhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_HyPol__resting.neomat'))
saFly_Lnop.flyID, saRest_Lnop.flyID
({0: 'fly1_G06xU29_20191104_cell1',
1: 'fly1_U29xG16_20191209_cell2',
2: 'fly3_X21_20200415_cell1',
3: 'fly3_x21_20200417_cell1',
4: 'fly2_X20_20210217_cell1'},
{0: 'fly1_G06xU29_20191104_cell1',
1: 'fly1_U29xG16_20191209_cell2',
2: 'fly3_X21_20200415_cell1',
3: 'fly3_x21_20200417_cell1',
4: 'fly2_X20_20210217_cell1'})
# saccading trials
saSacc_Lnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight',
'saccLeft_NoPol__CWTdetection_saccading.neomat'))
saSacc_Lnop.flyID
{0: 'fly1_G06xU29_20191104_cell1',
1: 'fly1_U29xG16_20191209_cell2',
2: 'fly3_X21_20200415_cell1',
3: 'fly3_x21_20200417_cell1',
4: 'fly2_X20_20210217_cell1'}
saRest_Lnop.average_flywise()
saSacc_Lnop.average_flywise()
saFly_Lnop.average_flywise()
# prettier, clean plotting, membrane potential
stimID = 3
# assert that we have equal datasets
restSet = set([v.lower() for _,v in saRest_Lnop.flyID.items()])
flySet = set([v.lower() for _,v in saSacc_Lnop.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Lnop.preStim+saSacc_Lnop.stimulusPhases_inframes[1]*saSacc_Lnop.meanFramePeriod)*saSacc_Lnop.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Lnop.samplingRate), int(maxLoom+postSacc*saSacc_Lnop.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Lnop.samplingRate-preSacc#-saSacc_Lnop[i].preStim
# d = saRest_Lnop.get_avg_data('ephys')
# s = saRest_Lnop.get_std_data('ephys')
d = saFly_Lnop.get_avg_data('ephys')
s = saFly_Lnop.get_std_data('ephys')
ax.plot(timeSlice, d[6][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[6]+s[6])[timeLims[0]:timeLims[1]], 100.*(d[6]-s[6])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Lnop.get_avg_data('ephys')
s = saRest_Lnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Lnop.get_avg_data('ephys')
s = saSacc_Lnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to left and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(['visual motion only','rest', 'saccades'])
ax.set_ylim(PLOT_Y_LIMS)
ax.set_xlim(PLOT_TIMEWINDOW)
#f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
saFly_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__flying.neomat'))
saRest_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__resting.neomat'))
saSacc_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
# saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight', 'newSorting',
# 'loomLeft_hyPol__CWTdetection_straight.neomat'))
saSacc_Rhyp.average_flywise()
saFly_Rhyp.average_flywise()
saRest_Rhyp.average_flywise()
# prettier, clean plotting, membrane potential
stimID = 0
# assert that we have equal datasets
restSet = set([v.lower() for k,v in saRest_Rhyp.flyID.items()])
flySet = set([v.lower() for k,v in saSacc_Rhyp.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rhyp.preStim+saSacc_Rhyp.stimulusPhases_inframes[1]*saSacc_Rhyp.meanFramePeriod)*saSacc_Rhyp.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rhyp.samplingRate), int(maxLoom+postSacc*saSacc_Rhyp.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rhyp.samplingRate-preSacc#-saSacc_Rhyp[i].preStim
d = saFly_Rhyp.get_avg_data('ephys')
s = saFly_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[6][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[6]+s[6])[timeLims[0]:timeLims[1]], 100.*(d[6]-s[6])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Rhyp.get_avg_data('ephys')
s = saRest_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Rhyp.get_avg_data('ephys')
s = saSacc_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS)
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
saFly = saFly_Rhyp
saSacc =
saRest =
saStrg = saStrg_Rhyp
# prettier, clean plotting, membrane potential
stimID = 0
# assert that we have equal datasets
restSet = set([v.lower() for k,v in saRest_Rhyp.flyID.items()])
flySet = set([v.lower() for k,v in saSacc_Rhyp.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.8 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rhyp.preStim+saSacc_Rhyp.stimulusPhases_inframes[1]*saSacc_Rhyp.meanFramePeriod)*saSacc_Rhyp.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rhyp.samplingRate), int(maxLoom+postSacc*saSacc_Rhyp.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rhyp.samplingRate-preSacc#-saSacc_Rhyp[i].preStim
d = saFly_Rhyp.get_avg_data('ephys')
s = saFly_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[6][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[6]+s[6])[timeLims[0]:timeLims[1]], 100.*(d[6]-s[6])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Rhyp.get_avg_data('ephys')
s = saRest_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Rhyp.get_avg_data('ephys')
s = saSacc_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS)
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
saFly_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__flying.neomat'))
saRest_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__resting.neomat'))
saSacc_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_saccading.neomat'))
# saStrg_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight', 'newSorting',
# 'loomLeft_noPol__CWTdetection_straight.neomat'))
saFly_Rnop.flyID, saRest_Rnop.flyID
({0: 'fly3_U29xG16_20200129_cell1',
1: 'fly1_X20a_20200316_cell1',
2: 'fly1_WTCS_20191106_cell1',
3: 'fly1_U29xG16_20200129_cell1',
4: 'fly1_G06xU29_20191104_cell1',
5: 'fly1_U29xG16_20191209_cell2',
6: 'fly3_X21_20200415_cell1',
7: 'fly3_x21_20200417_cell1',
8: 'fly2_X20_20210217_cell1'},
{0: 'fly3_U29xG16_20200129_cell1',
1: 'fly1_X20a_20200316_cell1',
2: 'fly1_WTCS_20191106_cell1',
3: 'fly1_U29xG16_20200129_cell1',
4: 'fly1_G06xU29_20191104_cell1',
5: 'fly1_U29xG16_20191209_cell2',
6: 'fly3_X21_20200415_cell1',
7: 'fly3_x21_20200417_cell1',
8: 'fly2_X20_20210217_cell1'})
saStrg_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_straight.neomat'))
saStrg_Rnop.clean_flyIDs()
saStrg_Rnop.reallocate_fly_indices()
saStrg_Rnop.write2neomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_straight.neomat'))
saSacc_Rnop.average_flywise()
saFly_Rnop.average_flywise()
saRest_Rnop.average_flywise()
# prettier, clean plotting, membrane potential
stimID = 0
# assert that we have equal datasets
restSet = set([v.lower() for k,v in saRest_Rnop.flyID.items()])
flySet = set([v.lower() for k,v in saSacc_Rnop.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rnop.preStim+saSacc_Rnop.stimulusPhases_inframes[1]*saSacc_Rnop.meanFramePeriod)*saSacc_Rnop.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rnop.samplingRate), int(maxLoom+postSacc*saSacc_Rnop.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rnop.samplingRate-preSacc#-saSacc_Rnop[i].preStim
d = saFly_Rnop.get_avg_data('ephys')
s = saFly_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[6][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[6]+s[6])[timeLims[0]:timeLims[1]], 100.*(d[6]-s[6])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Rnop.get_avg_data('ephys')
s = saRest_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Rnop.get_avg_data('ephys')
s = saSacc_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS)
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
Uses the same data and same code as Fig 4C, only difference is the stimID
saFly_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__flying.neomat'))
saRest_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__resting.neomat'))
saSacc_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
# saStrg_Rhyp = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight', 'newSorting',
# 'loomLeft_hyPol__CWTdetection_straight.neomat'))
saSacc_Rhyp.average_flywise()
saFly_Rhyp.average_flywise()
saRest_Rhyp.average_flywise()
# prettier, clean plotting, membrane potential
stimID = 1
# assert that we have equal datasets
restSet = set([v.lower() for k,v in saRest_Rhyp.flyID.items()])
flySet = set([v.lower() for k,v in saSacc_Rhyp.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rhyp.preStim+saSacc_Rhyp.stimulusPhases_inframes[1]*saSacc_Rhyp.meanFramePeriod)*saSacc_Rhyp.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rhyp.samplingRate), int(maxLoom+postSacc*saSacc_Rhyp.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rhyp.samplingRate-preSacc#-saSacc_Rhyp[i].preStim
d = saFly_Rhyp.get_avg_data('ephys')
s = saFly_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[7][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[7]+s[7])[timeLims[0]:timeLims[1]], 100.*(d[7]-s[7])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Rhyp.get_avg_data('ephys')
s = saRest_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Rhyp.get_avg_data('ephys')
s = saSacc_Rhyp.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS-3.8)# needs adjustment to fit, not changing interval size
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
saFly_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__flying.neomat'))
saRest_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__resting.neomat'))
saSacc_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_noPol__CWTdetection_saccading.neomat'))
# saStrg_Rnop = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
# 'Collections_saccade_straight', 'newSorting',
# 'loomLeft_noPol__CWTdetection_straight.neomat'))
saSacc_Rnop.average_flywise()
saFly_Rnop.average_flywise()
saRest_Rnop.average_flywise()
# prettier, clean plotting, membrane potential
stimID = 1
# assert that we have equal datasets
restSet = set([v.lower() for k,v in saRest_Rnop.flyID.items()])
flySet = set([v.lower() for k,v in saSacc_Rnop.flyID.items()])
assert(flySet == restSet)
f,ax=plt.subplots(1,1)
preSacc = 0.65 #s
postSacc = 0.5 #s
maxLoom = (saSacc_Rnop.preStim+saSacc_Rnop.stimulusPhases_inframes[1]*saSacc_Rnop.meanFramePeriod)*saSacc_Rnop.samplingRate
timeLims = [int(maxLoom-preSacc*saSacc_Rnop.samplingRate), int(maxLoom+postSacc*saSacc_Rnop.samplingRate)]
timeSlice = np.arange(timeLims[1]-timeLims[0])/saSacc_Rnop.samplingRate-preSacc#-saSacc_Rnop[i].preStim
d = saFly_Rnop.get_avg_data('ephys')
s = saFly_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[7][timeLims[0]:timeLims[1]]*100., color=optomotorColor)
ax.fill_between(timeSlice, 100.*(d[7]+s[7])[timeLims[0]:timeLims[1]], 100.*(d[7]-s[7])[timeLims[0]:timeLims[1]],
alpha=optomotorAlpha, color=optomotorColor)
d = saRest_Rnop.get_avg_data('ephys')
s = saRest_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[1])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[1], alpha=plotAlphas[1])
d = saSacc_Rnop.get_avg_data('ephys')
s = saSacc_Rnop.get_std_data('ephys')
ax.plot(timeSlice, d[stimID][timeLims[0]:timeLims[1]]*100., color=plotColors[2])
ax.fill_between(timeSlice, 100.*(d[stimID]+s[stimID])[timeLims[0]:timeLims[1]], 100.*(d[stimID]-s[stimID])[timeLims[0]:timeLims[1]],
color=plotColors[2], alpha=plotAlphas[2])
#ax.set_title('saccades to right and rotation right')
ax.set_ylabel('U [mV]')
ax.set_xlabel('time [s]')
ax.axvline(0, linestyle='--', color='k')
ax.legend(PLOT_LEGEND)
ax.set_ylim(PLOT_Y_LIMS-3.8)# needs adjustment to fit, not changing interval size
ax.set_xlim(PLOT_TIMEWINDOW)
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False, sizey=2.0, labely='2 mV', barwidth=3)
plt.tight_layout()
anlyz = lwb.Looming_with_background_Analyzer.retrieveFromNeomatfile(os.path.join(DATAPATH,
'Collections_saccade_straight', 'newSorting',
'loomLeft_hyPol__CWTdetection_saccading.neomat'))
# keep only one fly
rm = np.where(anlyz.fly_indices!=3)[0]
anlyz.remove_duds_byList(rm)
anlyz.clean_flyIDs()
anlyz.reallocate_fly_indices()
print(anlyz.flyID)
{0: 'fly1_WTCS_20190726_cell1'}
#choose a trial
trialID = 0 # they are actually all pretty nice
Plot raw L-R WSA trace and the derived PWC-filtered version of the signal.
yfit, xsteps, ysteps = helperLibrary.straightLineFiltering(
anlyz.segments[trialID].data[anlyz.channelmap['wingDiff']], 0.020*anlyz.samplingRate, postEnhancement=True )
f,ax=plt.subplots(1,1)
t = anlyz.timeax
d = anlyz.segments[trialID].data[anlyz.channelmap['wingDiff']]
ax.plot(t[:d.size], d, 'k-')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_ylabel('L-R WBA [°]')#
ax.set_xlabel('time from trial onset [s]')
ax.set_xticks([0.,0.5,1.0,1.5])
axins = ax.inset_axes([0.1, 0.55, 0.44, 0.5])
insetExtent = np.nonzero((t>0.44) & (t<0.68))[0]
axins.plot(t[insetExtent], d[insetExtent], 'k-', alpha=1.)
axins.set_ylim([0.185,0.26])
axins.set_xlim(t[insetExtent[[0,-1]]])
axins.set_yticks([0.2,0.25])
axins.set_xticks([0.45,0.65])
ax.indicate_inset_zoom(axins, edgecolor="black")
ax.plot(t[:yfit.size],yfit, '-', color=[31/255,119/255,180/255], alpha=1.)
axins.plot(t[insetExtent], yfit[insetExtent], '-', color=[31/255,119/255,180/255], alpha=1)
[<matplotlib.lines.Line2D at 0x18786ac2688>]
def central_diff_derivative(y, halfwidth=0.070, samplPeriod=0.020):
"""y is data to derivate, assumed to be np.array
halfwidth and samplPeriod shall use the same unit (e.g. seconds)"""
stepBack = int(halfwidth/samplPeriod+0.5)
stepFwd = int(halfwidth/samplPeriod)
return ( y[(stepFwd+stepBack):]-y[:-(stepBack+stepFwd)] )/((stepFwd+stepBack)*samplPeriod)
f,ax=plt.subplots(1,1)
ax2=ax.twinx()
medianFilterWidth = 3
diff_half_width = 0.05
detectSignal = signal.medfilt(ysteps,medianFilterWidth)
detectSignal = central_diff_derivative(detectSignal, diff_half_width, 0.02) #20 ms is Kinefly sampling period
ax.plot(ysteps, color='tab:blue')
ax2.plot(np.arange(int(diff_half_width/anlyz.kineflyPeriod+0.5),len(ysteps)-int(diff_half_width/anlyz.kineflyPeriod)) ,
detectSignal, color='tab:orange', linewidth=2)
# add picture of wavelet
ax.axhline(0,linestyle=':', color='k')
scales = 2.**np.arange(2, 3., 1./16.)
offset = 52
ax.plot(np.arange(-32,32.1)+offset, signal.ricker(65,scales[0]), color=[0,0.5,0], linestyle='-.', linewidth=2)
ax.plot(np.arange(-32,32.1)+offset, signal.ricker(65,scales[-1]), color=[0.6,0.5,0], linestyle='--',linewidth=2)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_ylabel('fitted L-R WBA [°]', color='tab:blue')
ax.tick_params(axis='y',labelcolor='tab:blue')
ax2.set_ylabel('$\delta_h\ [L-R\ WBA]$ [°/frame]', color='tab:orange')
ax2.tick_params(axis='y',labelcolor='tab:orange')
ax.set_xlabel('detection frame')
ax.spines['top'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
plt.tight_layout()
f,ax=plt.subplots(2,1)
scales = 2.**np.arange(2, 3., 1./16.)
detectSignal = signal.medfilt(ysteps,medianFilterWidth)
detectSignal = central_diff_derivative(detectSignal, diff_half_width, 0.02,)
# cwtSig = signal.cwt(ysteps,signal.ricker, scales)
cwtSig = signal.cwt(detectSignal,signal.ricker, scales)
# cwtSig,sj,freqs,coi,fft,fftfreqs = pycwt.cwt(detectSignal, 1, 1/16, s0=2,J=15, wavelet='mexicanhat')
# cwtSig= cwtSig.real
# ax[0].pcolormesh(cwtSig, cmap=plt.cm.coolwarm)
# im0=ax[0].pcolormesh(cwtSig, cmap=plt.cm.cividis, linewidth=0, rasterized=True)
# im0.set_edgecolor('face')
im0=ax[0].imshow(cwtSig, cmap=plt.cm.cividis, origin='lower')
cwtSig = signal.cwt(-detectSignal,signal.ricker, scales)
# im1=ax[1].pcolormesh(cwtSig, cmap=plt.cm.cividis, linewidth=0, rasterized=True)
# im1.set_edgecolor('face')
im1=ax[1].imshow(cwtSig, cmap=plt.cm.cividis, origin='lower')
# im0.set_clim([-0.5,0.5])
# im1.set_clim([-0.5,0.5])
# ax[0].set_yticklabels(['$2^{32/16}$','$2^{38/16}$','$2^{44/16}$', '$2^{50/16}$'])
# ax[0].set_yticklabels(['$2^2$','$2^{2+5/16}$','$2^{2+10/16}$', '$2^{2+15/16}$'])
# ax[0].set_yticks([0,8,16])
# ax[0].set_yticklabels(['$2^2$','$2^{2+8/16}$','$2^{2+16/16}$'])
ax[0].set_yticks([0,8,15])
ax[0].set_yticklabels(['$2^2$','$2^{2+8/16}$','$2^{2+15/16}$'])
# ax[1].set_yticks([0,8,16])
# ax[1].set_yticklabels(['$2^2$','$2^{2+8/16}$','$2^{2+16/16}$'])
ax[1].set_yticks([0,8,15])
ax[1].set_yticklabels(['$2^2$','$2^{2+8/16}$','$2^{2+15/16}$'])
# ax[0].set_title('CWT of signal')
ax[0].set_title('CWT of derivative')
# ax[1].set_title('CWT of inverse signal')
ax[1].set_title('CWT of negative derivative')
ax[1].set_xlabel('detection frame No')
# try 3d plot
# f=plt.figure()
# a0=f.add_subplot(projection='3d')
# a1=f.add_subplot(projection='3d')
# ax=[a0,a1]
# X,Y = np.meshgrid(np.arange(len(ysteps)), scales )
# a0.plot_surface(X,Y,cwtSig, cmap=plt.cm.coolwarm)
plt.tight_layout()
Needs additional function declarations to redo the saccade detection as it happens in the analysis algorithm. Disentagle overlapping point markers later in Inkscape.
def lowPass(y, samplingRate):
fThresh = 6.#Hz
#b,a=signal.butter(5,6.,fs=samplingRate, btype='low')
#ylow = signal.lfilter(b,a,y)
b,a=signal.butter(5,6.*1.5,fs=samplingRate, btype='low')
ylow = signal.filtfilt(b,a,y)
return ylow
def central_diff_derivative(y, halfwidth=0.070, samplPeriod=0.020):
"""y is data to derivate, assumed to be np.array
halfwidth and samplPeriod shall use the same unit (e.g. seconds)"""
stepBack = int(halfwidth/samplPeriod+0.5)
stepFwd = int(halfwidth/samplPeriod)
return ( y[(stepFwd+stepBack):]-y[:-(stepBack+stepFwd)] )/((stepFwd+stepBack)*samplPeriod)
def highpass_data(y, samplingRate):
# highpass data
b,a = signal.butter(5, 1., fs=samplingRate, btype='high')
yhigh = signal.lfilter(b,a,y)
return yhigh
def plotHighLowPass(y, yhigh, ylow=None):
f,ax=plt.subplots(1,1)
ax.plot(y)
ax.plot(yhigh)
if not ylow is None:
ax.plot(ylow)
ax.legend(['original','highpass', 'lowpass'])
return ax
def extractFlyStatistics(y, samplingRate, plot=True, **kwargs):
yhigh = highpass_data(y, samplingRate)
pdf, bins = np.histogram(yhigh, bins=50)
pdf = pdf/pdf.sum()
cdf = np.cumsum(pdf)
thresh = kwargs.get('thresh',0.95)
binCenters = (bins[:-1]+bins[1:])/2
minIdx = np.argmin(np.abs(cdf-thresh))
theta = binCenters[minIdx]
if plot:
f,ax=plt.subplots(1,1)
ax.plot(binCenters, pdf)
ax.plot(binCenters, cdf)
ax.plot(binCenters[minIdx],cdf[minIdx], 'go')
ax.plot(binCenters[minIdx],pdf[minIdx], 'ro')
ax.axhline(0.9, color='k', linestyle=':')
ax.axvline(theta, color='r')
return theta
def detectPeaks(theta, ydiff):
detectTheta = theta*6. #as stated in paper
peaksFound = signal.find_peaks(np.abs(ydiff), height=detectTheta)
return peaksFound
def plotDetection(ax,ydiff, detectTheta):
ax.plot(ydiff)
ax.axhline(detectTheta, color='k')
ax.axhline(0, color='g')
return
# find peak times
def findPeakTimes(ydiff, peaksFound):
## for each peak, move FORWARD till encountering closest to 0
peakTime = np.full_like(peaksFound, np.nan, dtype=int)
for j,p in enumerate(peaksFound):
#compare the signs, because peaks can be + or - (at 0 it switches)
sameSign = True
i=0
while sameSign>0:
i+=1
if p+i >=ydiff.size:
i -=1 #plug right value for line after while loop
break
sameSign = (ydiff[p+i]*ydiff[p+i-1])>0.
if p+i <ydiff.size:
thisPeak = p+i-1 + np.argmin(np.abs(ydiff[[p+i-1,p+i]]))
else:
thisPeak = ydiff.size
peakTime[j] = int(thisPeak)
return peakTime
# find onset times
def findOnsetTime(ydiff, peaksFound):
## for each peak, move BACKWARD till encountering closest to 0
onsetTime = np.zeros_like(peaksFound, dtype=int)
for j,p in enumerate(peaksFound):
sameSign = True
i=-0
while sameSign>0:
i-=1
if p+i<0:
i+=1 #plug right value for line after while loop
break
sameSign = (ydiff[p+i]*ydiff[p+i+1])>0. #compare the signs, because peaks can be + or - (at 0 it switches)
if p+i >=0:
thisOnset = p+i + np.argmin(np.abs(ydiff[[p+i,p+i+1]]))
else:
thisOnset=0
onsetTime[j] = int(thisOnset)
return onsetTime
# correct onset times for error introduced by lowpass filtering
def correctOnsetTime(y,ylow, ydiff, peaksFound, onsetTime):
actualOnset = np.full_like(onsetTime, np.nan, dtype=int)
for j,o in enumerate(onsetTime):
## use value of low-passed L-R WBA at putative onset as threshold
thresh = ylow[o]
## look forward in original signal from putative onset and find last sample before threshold
ytemp = y-thresh #this turns it into finding a 0-crossing (like above)
### shortcut to deal with cases where actual onset is not included in trace
p = peaksFound[j]
if not any(np.diff(np.sign(ytemp[o:p]))):
#no change in sign means no threshold crossing
actualOnset[j] = o #keep putataive one
continue
sameSign=True
i=0
while sameSign:
i+=1
sameSign = (ytemp[o+i]*ytemp[o+i-1])>0.
## define this as actual onset
actualOnset[j] = o+1-1
return actualOnset
## get saccade amplitudes
def getSaccadeAmplitude(y, peaksFound, actualOnset, samplingRate):
#subtracted the mean L-R WBA signal 50-ms baseline interval before saccade onset
#from mean L-R WBA signal in +-15-ms windows surrounding the peak
saccAmp = np.full(peaksFound.size, np.nan, dtype=np.float)
for i in range(len(actualOnset)):
baseline = np.mean(y[actualOnset[i]-int(0.06*samplingRate) : actualOnset[i]])
if not np.isnan(baseline):
saccAmp[i] = np.mean( y[peaksFound[i]-int(0.03*samplingRate) : peaksFound[i]+int(0.03*samplingRate)] ) - baseline
return saccAmp
def plotPeaks(y, actualOnset, peakTime, saccAmp, ax=None):
if ax is None:
f,ax=plt.subplots(1,1)
ax.plot(y)
ax.plot(actualOnset, y[actualOnset], 'r>')
ax.plot(peakTime, y[peakTime], 'k^')
ax.plot(peakTime, saccAmp, 'ko')
return
# eliminate putative saccade if maximal deviation before onset is>55% of amplitude
def eliminateByDeviation(y, actualOnset, saccAmp, samplingRate, **kwargs):
window = kwargs.get('window',0.201)
deviationLimit = kwargs.get('deviationLimit',0.55)
toEliminate = []
for i in range(actualOnset.size):
# preselect to make sure the indexing is feasible
chunk = y[actualOnset[i]-int(window*samplingRate) : actualOnset[i]]
if actualOnset[i]-int(window*samplingRate) <0:
toEliminate.append(i)
continue
deviation = chunk.max()-chunk.min()
if deviation>deviationLimit*abs(saccAmp[i]):
toEliminate.append(i)
return toEliminate
def eliminateByPeakProximity(onset,saccAmp,*, samplingRate=10000, timeDist:float=0.5) ->list:
'''eliminate peaks that are too close to one another.
This is intended particularly for when they go in opposite direction,
such that the second could well be end of the first
Note: we don't need absolute times because operating inside segments
This is instance three, considering only distances between saccades if they are in opposite directions'''
#make sure they are sorted
onset.sort()
peakTime.sort()
# iterate forward, check distance to previous
toEliminate = []
for i in range(1,onset.size):
if saccAmp[i]*saccAmp[i-1] >0:#sign check, same saccade direction
continue
current=onset[i]
prev = onset[i-1]
# if it's too short, eliminate
assert(not current-prev<0)
if current-prev < timeDist*samplingRate :
toEliminate.append(i)
return toEliminate
def eliminateBySaccAmplitude(saccAmp, thresh)->list:
toEliminate=[]
for i in range(0,saccAmp.size):
if abs(saccAmp[i])< thresh:
toEliminate.append(i)
return toEliminate
def eliminatePutativePersistent(cwtpeaks,peakTime,maxPeakDelay,onset,maxOnsetDelay, samplingPeriod):
# eliminate if the onset or peak are too far away from the peak in derivative
toEliminate = []
for i in range(cwtpeaks.size):
if cwtpeaks[i]-onset[i]>maxOnsetDelay/samplingPeriod:#don't know yet how big this hsould be
toEliminate.append(i)
elif peakTime[i]-cwtpeaks[i]>maxPeakDelay/samplingPeriod:
toEliminate.append(i)
return toEliminate
def align_yaxis(ax1, v1, ax2, v2):
"""adjust ax2 ylimit so that v2 in ax2 is aligned to v1 in ax1"""
_, y1 = ax1.transData.transform((0, v1))
_, y2 = ax2.transData.transform((0, v2))
inv = ax2.transData.inverted()
_, dy = inv.transform((0, 0)) - inv.transform((0, y1-y2))
miny, maxy = ax2.get_ylim()
ax2.set_ylim(miny+dy, maxy+dy)
return
#saccDataCWT = {}
scales = 2.**np.arange(2, 3., 1./16.)
diff_half_width = 0.05 #0.05
medianFilterWidth = 3
cwtpeaks_window_size = 20
deviationWindow = 0.16 #0.16
deviationLimit = 0.65 # was 0.55 for spont
saccAmpThreshold = 15.*np.pi/180
minPeakInterval = 0.1
cascadeFiltering = True
maxPeakDelay = 1.2#s
maxOnsetDelay = 1.1#s
f,ax=plt.subplots(1,1)
ax2 = ax.twinx()
for i in [trialID]: #replace iteration by length-1-list for this plot
_,_,ysteps = helperLibrary.straightLineFiltering(anlyz.segments[trialID].data[anlyz.channelmap['wingDiff']],
0.020*anlyz.samplingRate, postEnhancement=True )
y=ysteps
detectSignal = signal.medfilt(y,medianFilterWidth)
detectSignal = central_diff_derivative(detectSignal, diff_half_width, 0.02)
ax2.plot(range(1,y.size), np.diff(signal.medfilt(y,medianFilterWidth)), color='tab:green',
label='differential', linestyle='--', linewidth=2.5)
ax2.axhline(0, color='k',linestyle=':')
peaks = signal.find_peaks_cwt(detectSignal, scales, gap_thresh=1,
window_size=cwtpeaks_window_size, noise_perc=5)
peaks2 = signal.find_peaks_cwt(-detectSignal, scales, gap_thresh=1,
window_size=cwtpeaks_window_size, noise_perc=5)
peaks = np.concatenate((peaks,peaks2))
peaks.sort() #sorts in-place
peaks = peaks.astype('int')
peaks += int(diff_half_width/anlyz.kineflyPeriod)
peakTime = findPeakTimes(np.diff(signal.medfilt(y,medianFilterWidth)), peaks-1) +1 #plus 1 to correct diff shift
onsetTime = findOnsetTime(np.diff(signal.medfilt(y,medianFilterWidth)), peaks-1) +1
saccAmp = getSaccadeAmplitude(y, peakTime, onsetTime, 1/anlyz.kineflyPeriod)
# plot detection before filtering
ax.plot(y, linewidth=2)
ax.plot(peaks, y[peaks], 'x',color='tab:orange', markersize=10, markeredgewidth=2, label='differential peak')
ax.plot(peakTime, y[peakTime], '^', color='k', markersize=8, label='candidate peak')
ax.plot(onsetTime, y[onsetTime], '>', color='k', markersize=8, label='candidate onset')
#start filtering stage
keep = np.ones_like(peaks,dtype=bool)
elim = eliminateBySaccAmplitude(saccAmp, saccAmpThreshold)
keep[elim]=False
if cascadeFiltering:
peaks = peaks[keep]
peakTime = peakTime[keep]
onsetTime = onsetTime[keep]
saccAmp = saccAmp[keep]
keep = np.ones_like(peaks,dtype=bool)
elim = eliminateByPeakProximity(onsetTime, saccAmp, samplingRate=1/anlyz.kineflyPeriod,
timeDist=minPeakInterval)
keep[elim]=False
#combine with next step
#these two steps are similar, therefore combine together
elim = eliminateByDeviation(y, onsetTime, saccAmp, 1/anlyz.kineflyPeriod,
window=deviationWindow, deviationLimit=deviationLimit)
keep[elim]=False
if cascadeFiltering:
peaks = peaks[keep]
peakTime = peakTime[keep]
onsetTime = onsetTime[keep]
saccAmp = saccAmp[keep]
keep = np.ones_like(peaks,dtype=bool)
elim = eliminatePutativePersistent(peaks, peakTime,maxPeakDelay, onsetTime,maxOnsetDelay, anlyz.kineflyPeriod)
keep[elim]=False
ax.plot(peaks[keep], y[peaks[keep]], 'o',color='r', markerfacecolor="None",
markeredgewidth=2, markersize=11)
ax.plot(peakTime[keep], y[peakTime[keep]], '^',color='r',markeredgewidth=2,
markerfacecolor="None", markersize=9, label='definitive peak')
ax.plot(onsetTime[keep], y[onsetTime[keep]], '>',color='r',markeredgewidth=2,
markerfacecolor="None", markersize=9, label='definitive onset')
ax.set_title(trialID)
peaks = peaks[keep]
peakTime = peakTime[keep]
onsetTime = onsetTime[keep]
saccAmp = saccAmp[keep]
# saccDataCWT[i] = dict([('peakTime',peakTime),('onset',onsetTime), ('saccAmp',saccAmp), ('CWTpeaks',peaks)])
"""detectionParams = {'method':'CWT',
'scales': scales,
'diff_half_width': diff_half_width,
'medianFilterWidth': medianFilterWidth,
'CWTpeaks_window_size':cwtpeaks_window_size,
'deviationWindow':deviationWindow,
'deviationLimit': deviationLimit,
'saccAmpThreshold': saccAmpThreshold,
'minPeakInterval':minPeakInterval,
'maxPeakDelay': maxPeakDelay,
'maxOnsetDelay': maxOnsetDelay,
'cascadeFiltering': cascadeFiltering
}
"""
C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:131: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
"detectionParams = {'method':'CWT',\n 'scales': scales,\n 'diff_half_width': diff_half_width,\n 'medianFilterWidth': medianFilterWidth,\n 'CWTpeaks_window_size':cwtpeaks_window_size,\n 'deviationWindow':deviationWindow,\n 'deviationLimit': deviationLimit,\n 'saccAmpThreshold': saccAmpThreshold,\n 'minPeakInterval':minPeakInterval,\n 'maxPeakDelay': maxPeakDelay,\n 'maxOnsetDelay': maxOnsetDelay,\n 'cascadeFiltering': cascadeFiltering\n }\n"
flist = glob.glob(os.path.join(DATAPATH,'spontaneousSaccades','oldData','CWTdetection','*.neomat'))
flist
['D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20111222_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120104_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120104_cell2_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120105_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120105_cell2_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120207_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120217_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120418_cell2_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120501_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120504_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120508_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120524_cell2_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20120605_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20130211_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20130424_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20130502_cell2_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20130506_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\oldData\\CWTdetection\\fly1_wtcs_20130513_cell1_oldData_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly4_x21_20200420_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly3_x21_20200415_cell1_looming_wb_spontSaccades_cwtDetection.neomat', 'D:\\RecorderData\\spontaneousSaccades\\CWTdetection\\fly3_x21_20200417_cell1_looming_wb_spontSaccades_cwtDetection.neomat']
allSaccLeftMeanEphys = {}
allSaccLeftMeanWBA = {}
allSaccRightMeanEphys = {}
allSaccRightMeanWBA = {}
onset = 0.15
for fn in flist[:]:
da = DummyAnalyzer.retrieveFromNeomatfile(fn)
if 'oldData' not in fn: # conversion from radian to degree for data recorded in this study
print('multiplying WBA at',fn)
da.saccLeftMeanWBA *= (180/np.pi)
da.saccRightMeanWBA *= (180/np.pi)
else: # old data is from a different rig, and needs another conversion factor to obtain degrees
da.saccLeftMeanWBA /= 1.7
da.saccRightMeanWBA /= 1.7
if max(abs(da.saccLeftMeanEphys)) <0.1: #same idea, deal with different conversion factors
print('multiplying Vm at',fn)
da.saccLeftMeanEphys *= 100.
da.saccRightMeanEphys *= 100.
if abs(da.samplingRate -10000) >10:
raise Exception('unexpected sampling rate')
# gather part of data that matters for plots, i.e. 150ms before onset and the 500ms after
allSaccLeftMeanEphys[da.flyID[0]] = da.saccLeftMeanEphys[-int((0.5+onset)*da.samplingRate):]
allSaccRightMeanEphys[da.flyID[0]] = da.saccRightMeanEphys[-int((0.5+onset)*da.samplingRate):]
allSaccLeftMeanWBA[da.flyID[0]] = da.saccLeftMeanWBA[-int((0.5+onset)*da.samplingRate):]
allSaccRightMeanWBA[da.flyID[0]] = da.saccRightMeanWBA[-int((0.5+onset)*da.samplingRate):]
multiplying WBA at D:\RecorderData\spontaneousSaccades\CWTdetection\fly4_x21_20200420_cell1_looming_wb_spontSaccades_cwtDetection.neomat multiplying Vm at D:\RecorderData\spontaneousSaccades\CWTdetection\fly4_x21_20200420_cell1_looming_wb_spontSaccades_cwtDetection.neomat multiplying WBA at D:\RecorderData\spontaneousSaccades\CWTdetection\fly3_x21_20200415_cell1_looming_wb_spontSaccades_cwtDetection.neomat multiplying WBA at D:\RecorderData\spontaneousSaccades\CWTdetection\fly3_x21_20200417_cell1_looming_wb_spontSaccades_cwtDetection.neomat multiplying Vm at D:\RecorderData\spontaneousSaccades\CWTdetection\fly3_x21_20200417_cell1_looming_wb_spontSaccades_cwtDetection.neomat
subtype_dict = {'HSN': ('fly1_wtcs_20111222_cell1',
'fly1_wtcs_20120207_cell1',
'fly1_wtcs_20120504_cell1',
'fly1_wtcs_20120508_cell1',
'fly1_wtcs_20120418_cell2',
'fly1_wtcs_20130506_cell1',
),
'HSE': ('fly1_wtcs_20120104_cell1',
'fly1_wtcs_20120105_cell1',
'fly1_wtcs_20120501_cell1',
'fly1_wtcs_20120105_cell2',
'fly1_wtcs_20120524_cell2',
'fly1_wtcs_20130502_cell2',
'fly1_wtcs_20130513_cell1'),
'HSS': ('fly1_wtcs_20120605_cell1',
'fly1_wtcs_20120104_cell2',
'fly1_wtcs_20130211_cell1',
'fly1_wtcs_20130424_cell1',
)
}
f,ax=plt.subplots(1,1, figsize=(7,5))
# BEGIN WITH HSE
plotlist = subtype_dict['HSE']
print('HSE')
plotcolors = plt.cm.YlGn(np.linspace(0.2,0.7,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanWBA[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanWBA[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanWBA[fid][:minLen] for fid in plotlist])
avg_HSE = np.mean(alldata, axis=0)
std_HSE = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSE.size)/10000.-0.15,avg_HSE, color=plt.cm.YlGn(0.85), linewidth=2)
# THEN HSS
plotlist = subtype_dict['HSS']
print('HSS')
plotcolors = plt.cm.pink(np.linspace(0.3,0.4,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanWBA[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanWBA[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanWBA[fid][:minLen] for fid in plotlist])
avg_HSS = np.mean(alldata, axis=0)
std_HSS = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSS.size)/10000.-0.15,avg_HSS, color=plt.cm.Purples(0.98), linewidth=2)
# LAST COME HSN
plotlist = subtype_dict['HSN']
print('HSN')
plotcolors = plt.cm.Purples(np.linspace(0.2,0.7,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanWBA[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanWBA[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanWBA[fid][:minLen] for fid in plotlist])
avg_HSN = np.mean(alldata, axis=0)
std_HSN = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSN.size)/10000.-0.15,avg_HSN, color=plt.cm.Purples(0.98), linewidth=2)
# FINALLY ALL AVERAGE LINES
ax.plot(np.arange(avg_HSS.size)/10000.-0.15,avg_HSS, color=plt.cm.pink(0.2), linewidth=3, label='mean HSS')
ax.plot(np.arange(avg_HSE.size)/10000.-0.15,avg_HSE, color=plt.cm.YlGn(0.85), linewidth=3, label='mean HSE')
ax.plot(np.arange(avg_HSN.size)/10000.-0.15,avg_HSN, color=plt.cm.Purples(0.98), linewidth=3, label='mean HSN')
# ax.legend()
# sugar
ax.axvline(0, linestyle='--', color='k')
ax.set_xlim([-0.15,0.48])
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False,
sizey=5.0, labely='5°', barwidth=3)
plt.tight_layout()
HSE fly1_wtcs_20120104_cell1 fly1_wtcs_20120105_cell1 fly1_wtcs_20120501_cell1 fly1_wtcs_20120105_cell2 fly1_wtcs_20120524_cell2 fly1_wtcs_20130502_cell2 fly1_wtcs_20130513_cell1 HSS fly1_wtcs_20120605_cell1 fly1_wtcs_20120104_cell2 fly1_wtcs_20130211_cell1 fly1_wtcs_20130424_cell1 fly3_x21_20200415_cell1 fly3_x21_20200417_cell1 HSN fly1_wtcs_20111222_cell1 fly1_wtcs_20120207_cell1 fly1_wtcs_20120504_cell1 fly1_wtcs_20120508_cell1 fly1_wtcs_20120418_cell2 fly1_wtcs_20130506_cell1 fly4_x21_20200420_cell1
This is essentially the same plot as in Fig S2B, only replacing the L-R WBA with ephys
f,ax=plt.subplots(1,1, figsize=(7,5))
# BEGIN WITH HSE
plotlist = subtype_dict['HSE']
print('HSE')
plotcolors = plt.cm.YlGn(np.linspace(0.2,0.7,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanEphys[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanEphys[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanEphys[fid][:minLen] for fid in plotlist])
avg_HSE = np.mean(alldata, axis=0)
std_HSE = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSE.size)/10000.-0.15,avg_HSE, color=plt.cm.YlGn(0.85), linewidth=2)
# THEN HSS
plotlist = subtype_dict['HSS']
print('HSS')
plotcolors = plt.cm.pink(np.linspace(0.3,0.4,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanEphys[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanEphys[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanEphys[fid][:minLen] for fid in plotlist])
avg_HSS = np.mean(alldata, axis=0)
std_HSS = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSS.size)/10000.-0.15,avg_HSS, color=plt.cm.Purples(0.98), linewidth=2)
# LAST HSN
plotlist = subtype_dict['HSN']
print('HSN')
plotcolors = plt.cm.Purples(np.linspace(0.2,0.7,len(plotlist)))
for i,k in enumerate(plotlist):
print(k)
data = allSaccLeftMeanEphys[k]
ax.plot(np.arange(data.size)/10000.-0.15, data, color=plotcolors[i], linewidth=1., label=k)
# add average line
minLen = min([allSaccLeftMeanEphys[fid].size for fid in plotlist])
alldata = np.array([allSaccLeftMeanEphys[fid][:minLen] for fid in plotlist])
avg_HSN = np.mean(alldata, axis=0)
std_HSN = np.std(alldata, axis=0)/len(plotlist)
# ax.plot(np.arange(avg_HSN.size)/10000.-0.15,avg_HSN, color=plt.cm.Purples(0.98), linewidth=2)
# ADD AVERAGE LINES (plotted here so they end up on top)
ax.plot(np.arange(avg_HSS.size)/10000.-0.15,avg_HSS, color=plt.cm.pink(0.2), linewidth=3, label='mean HSS')
ax.plot(np.arange(avg_HSE.size)/10000.-0.15,avg_HSE, color=plt.cm.YlGn(0.85), linewidth=3, label='mean HSE')
ax.plot(np.arange(avg_HSN.size)/10000.-0.15,avg_HSN, color=plt.cm.Purples(0.99), linewidth=3, label='mean HSN')
# ax.legend()
# some sugar
ax.axvline(0, linestyle='--', color='k')
ax.set_xlim([-0.15,0.48])
# f.set_size_inches(85/25.4, 50/25.4)
add_scalebar(ax, matchx=False, sizex=0.1, labelx='100 ms', matchy=False,
sizey=1.0, labely='1 mV', barwidth=3)
plt.tight_layout()
HSE fly1_wtcs_20120104_cell1 fly1_wtcs_20120105_cell1 fly1_wtcs_20120501_cell1 fly1_wtcs_20120105_cell2 fly1_wtcs_20120524_cell2 fly1_wtcs_20130502_cell2 fly1_wtcs_20130513_cell1 HSS fly1_wtcs_20120605_cell1 fly1_wtcs_20120104_cell2 fly1_wtcs_20130211_cell1 fly1_wtcs_20130424_cell1 fly3_x21_20200415_cell1 fly3_x21_20200417_cell1 HSN fly1_wtcs_20111222_cell1 fly1_wtcs_20120207_cell1 fly1_wtcs_20120504_cell1 fly1_wtcs_20120508_cell1 fly1_wtcs_20120418_cell2 fly1_wtcs_20130506_cell1 fly4_x21_20200420_cell1
# extract maximal L-R WBA and membrane potential amplitudes in a 200ms window after saccade onset put in dict structure
# refer to previosu cell to avoid mistakes
polAmp_dict={}
saccAmp_dict={}
onset=0.15 #150ms into datapiece, because of way we loaded data at beginning of Figure S2
for flyid in allSaccRightMeanEphys.keys():
d=allSaccLeftMeanEphys[flyid]
d = ndimage.gaussian_filter1d(d, 60)
ix = int(np.argmax(signal.medfilt(np.abs(d[int(onset*10000):int((onset+0.2)*10000)]), 101))+0.15*10000)
polAmp_dict[flyid] = d[ix]
d=allSaccLeftMeanWBA[flyid]
d = ndimage.gaussian_filter1d(d, 60)
ix = int(np.argmax(signal.medfilt(np.abs(d[int(onset*10000):int((onset+0.2)*10000)]), 101))+0.15*10000)
saccAmp_dict[flyid] = d[ix]
#statisticType = 'median-quartile'
statisticType = 'mean-sd'
f,ax=plt.subplots(1,1, figsize=(5,4))
#HSE
plotlist = subtype_dict['HSE']
# remove the problematic one(s)
# plotlist = [flyid for flyid in subtype_dict['HSE'] if not flyid in behavioral_outliers]
plotcolors = plt.cm.YlGn(np.linspace(0.2,0.7,len(plotlist)))
i=0
print('HSE')
for i,flyid in enumerate(plotlist):
print(flyid)
ax.scatter(saccAmp_dict[flyid], polAmp_dict[flyid], color=plotcolors[i], marker='o')
i+=1
# calculate average dot
if statisticType=='mean-sd':
avg_HSE_wba = np.mean( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSE_ephys = np.mean( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSE_wba = np.std( [saccAmp_dict[flyid] for flyid in plotlist])
dev_HSE_ephys = np.std( [polAmp_dict[flyid] for flyid in plotlist])
elif statisticType=='median-quartile':# median and quartile
avg_HSE_wba = np.median( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSE_ephys = np.median( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSE_wba = np.quantile( [saccAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
dev_HSE_ephys = np.quantile( [polAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
else:
raise Exception('unexpected key word for statisticType')
# HSS
plotlist = subtype_dict['HSS']
# plotlist = [flyid for flyid in subtype_dict['HSS'] if not flyid in behavioral_outliers]
plotcolors = plt.cm.pink(np.linspace(0.3,0.4,len(plotlist)))
i=0
print('HSS')
for i,flyid in enumerate(plotlist):
print(flyid)
ax.plot(saccAmp_dict[flyid], polAmp_dict[flyid], color=plotcolors[i], marker='o')
i+=1
if statisticType=='mean-sd':
avg_HSS_wba = np.mean( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSS_ephys = np.mean( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSS_wba = np.std( [saccAmp_dict[flyid] for flyid in plotlist])
dev_HSS_ephys = np.std( [polAmp_dict[flyid] for flyid in plotlist])
elif statisticType=='median-quartile':# median and quartile
avg_HSS_wba = np.median( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSS_ephys = np.median( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSS_wba = np.quantile( [saccAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
dev_HSS_ephys = np.quantile( [polAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
else:
raise Exception('unexpected key word for statisticType')
# HSN
plotlist = subtype_dict['HSN']
# plotlist = [flyid for flyid in subtype_dict['HSN'] if not flyid in behavioral_outliers]
plotcolors = plt.cm.Purples(np.linspace(0.2,0.7,len(plotlist)))
i=0
print('HSN')
for i,flyid in enumerate(plotlist):
print(flyid)
ax.plot(saccAmp_dict[flyid], polAmp_dict[flyid], color=plotcolors[i], marker='o')
i+=1
if statisticType=='mean-sd':
avg_HSN_wba = np.mean( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSN_ephys = np.mean( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSN_wba = np.std( [saccAmp_dict[flyid] for flyid in plotlist])
dev_HSN_ephys = np.std( [polAmp_dict[flyid] for flyid in plotlist])
elif statisticType=='median-quartile':# median and quartile
avg_HSN_wba = np.median( [saccAmp_dict[flyid] for flyid in plotlist])
avg_HSN_ephys = np.median( [polAmp_dict[flyid] for flyid in plotlist])
dev_HSN_wba = np.quantile( [saccAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
dev_HSN_ephys = np.quantile( [polAmp_dict[flyid] for flyid in plotlist], [0.25,0.75])
else:
raise Exception('unexpected key word for statisticType')
#sugar
ax.set_xlabel('mean saccade amplitude [°]')
ax.set_ylabel('mean saccade-associated potential [mV]')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set_xticks(range(-22,-9,4))
ax.set_yticks(range(-2,5,2))
plt.tight_layout()
# add average dots
ax.plot(avg_HSE_wba, avg_HSE_ephys, color=plt.cm.YlGn(0.85),
marker='s',markersize=10,markeredgecolor='k',linewidth=2)
ax.plot(avg_HSS_wba, avg_HSS_ephys, color=plt.cm.pink(0.2),
marker='s',markersize=10,markeredgecolor='k',linewidth=2)
ax.plot(avg_HSN_wba, avg_HSN_ephys, color=plt.cm.Purples(0.98),
marker='s',markersize=10,markeredgecolor='k',linewidth=2)
# plotting mean +- sd
ax.errorbar(avg_HSE_wba, avg_HSE_ephys, xerr=dev_HSE_wba, yerr=dev_HSE_ephys, color=plt.cm.YlGn(0.85),
marker='s',markersize=10,markeredgecolor='k',linewidth=3)
ax.errorbar(avg_HSS_wba, avg_HSS_ephys, xerr=dev_HSS_wba, yerr=dev_HSS_ephys, color=plt.cm.pink(0.2),
marker='s',markersize=10,markeredgecolor='k',linewidth=3)
ax.errorbar(avg_HSN_wba, avg_HSN_ephys, xerr=dev_HSN_wba, yerr=dev_HSN_ephys, color=plt.cm.Purples(0.98),
marker='s',markersize=10,markeredgecolor='k',linewidth=3)
HSE fly1_wtcs_20120104_cell1 fly1_wtcs_20120105_cell1 fly1_wtcs_20120501_cell1 fly1_wtcs_20120105_cell2 fly1_wtcs_20120524_cell2 fly1_wtcs_20130502_cell2 fly1_wtcs_20130513_cell1 HSS fly1_wtcs_20120605_cell1 fly1_wtcs_20120104_cell2 fly1_wtcs_20130211_cell1 fly1_wtcs_20130424_cell1 fly3_x21_20200415_cell1 fly3_x21_20200417_cell1 HSN fly1_wtcs_20111222_cell1 fly1_wtcs_20120207_cell1 fly1_wtcs_20120504_cell1 fly1_wtcs_20120508_cell1 fly1_wtcs_20120418_cell2 fly1_wtcs_20130506_cell1 fly4_x21_20200420_cell1
<ErrorbarContainer object of 3 artists>